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REAL-TIME  IMPLEMENTATION  OF  THE  KALMAN  FILTER 
FOR  TRAJECTORY  ESTIMATION 

I  INTRODUCTION 

The  purpose  of  this  technical  memorandum  is  to  describe  the  results 
of  a  study  concerned  with  real-time  implementation  of  the  Kalman  filter 
for  estimating  ballistic  trajectories.  As  shown  in  previous  studies1 »2* 
performed  by  the  Information  and  Control  Laboratory  in  support  of  the 
NIKE-X  system,  the  performance  of  the  Kalman  filter  for  the  trajectory 
estimation  task  is  excellent;  however,  the  computation  time  required 
by  this  filter  is  relatively  large.  On  the  Univac  1108  computer,  the  com¬ 
putation  time  per  filter  iteration  is  approximately  0.010  sec  for  endo- 
atmospheric  estimation;  hence,  with  a  data  rate  of  20  measurements  per 
second,  the  filter  runs  about  5  times  faster  than  real  time.  For  exo- 
atmospheric  estimation  (in  which  the  filter  has  6  states  instead  of  7, 
and  is  programmed  in  double  precision— see  Sec.  VI),  the  filter's  itera¬ 
tion  time  is  approximately  0.011  sec;  hence,  with  a  data  rate  of  one 
measurement  per  second,  the  filter  runs  about  100  times  faster  than  real 
time. 

For  real-time  estimation  of  multiple  targets,  the  computation  speed 
of  the  Kalman  filter  may  be  inadequate  (at  least  for  endoatmospheric 
trajectory  estimation).  Thus,  it  is  essential  to  improve  the  computa¬ 
tional  efficiency  of  the  filtering  algorithm;  this  should  be  done  with 
a  minimum  amount  of  degradation  in  filter  performance. 

In  Sec.  II,  the  problem  formulation  is  presented;  this  consists 
of  the  equations  of  motion  for  the  ballistic  reentry  vehicle  and  the 
radar  measurement  equations.  The  equations  for  the  extended  Kalman 
filter  and  a  flow  diagram  depicting  its  operation  are  given  in  Sec.  III. 

In  Secs.  Ill  and  IV  several  approaches  that  may  be  used  to  modify 
the  Kalman  filtering  algorithm  in  order  to  reduce  the  computational 
requirements  are  described;  these  include 
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(1)  Use  of  a  precomputed  approximation  to  the  weighting  matrix 
in  the  filter  equations  (this  also  eliminates  the  necessity 
of  calculating  the  covariance  matrices);  i.e.,  the  weighting 
matrix  is  precomputed  off-line,  rather  than  calculated  recur¬ 
sively  on-line. 

(2)  Increasing  the  sample  Interval  (data  processing  interval) 

At,  the  computation  speed  of  the  filter  relative  to  real 
time  increases  as  this  parameter  increases. 

(3)  Updating  the  weighting  matrix  and  the  covariance  matrices  at 

intervals  of  At  sec,  where  ~  =  N  2;  1,  i.e.,  the  weighting 

At 

matrix  is  only  updated  every  N  iterations  and  is  held  con¬ 
stant  in  the  filter  equations  for  these  iterations,  before 
being  updated  anew.  This  computational  algorithm  will  be  refer¬ 
red  to  as  the  piecewise-recursive  Kalman  filter. 

The  first  appx’oach  was  studied  extensively  in  Ref.  1;  the  other 
two  approaches  will  be  investigated  in  considerable  detail  in  this  report. 
These  approaches  were  evaluated  and  compared  by  estimating  the  states  of 
a  ballistic  trajectory  based  on  simulated  radar  measurements  that  are 
corrupted  by  noise;  the  simulated  trajectories  were  generated  from  an 
accurate  model  of  reentry  dynamics.  The  numerical  results  that  were 
obtained  in  this  study  are  presented  in  Sec.  V  (for  endoatmospheric 
estimation)  and  Sec.  VI  (for  exoatmospheric  estimation).  Estimation 
errors  in  position,  velocity,  and  acceleration,  and  the  actual  and  estimated 
g  are  plotted  as  functions  of  time  for  the  various  computer  runs.  In 
addition,  the  computation  times  (on  the  Univac  1108)  per  iteration  for 
the  different  filter  algorithms  are  presented. 

For  real-time  estimation  of  multiple  targets  ,  the  piecewise-recur¬ 
sive  Kalman  filter  is  the  most  promising  approach  of  those  considered  in 
this  report.  For  endoatmospheric  estimation,  it  is  possible  to  filter 
measurements  at  a  computational  speed  (on  the  Univac  1108)  that  is  approxi¬ 
mately  20  or  25  times  faster  than  real  time  (for  At  =  0.05  sec,  At  =  0.5 
sec  or  At  =  0.10  sec,  At  =  0.5  sec,  respectively)  and  yet  obtain  accuracy 
approaching  that  of  the  fully  implemented  Kalman  filter;  for  exoatmospheric 
estimation,  the  computational  speed  (on  the  Univac  1108)  is  approximately 
500  times  faster  than  real  time  (for  At  =  1.0  sec,  At  =  20  sec). 
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The  last  section  of  the  report  summarizes  the  results  of  this 
study  and  indicates  areas  for  future  work. 


II  PROBLEM  FORMULATION 


A.  Equations  of  Motion  for  the  Ballistic  Reentry  Vehicle 

A  model  for  the  equations  of  motion  for  the  ballistic  reentry 
vehicle  is  taken  to  be* 

2  1  P&o  x 

x  =  -  2<i)(z  cos  t  -  y  sin  {,)  +  u>  x  -  —  Vx  -  g  — - 

2  P  T1 

2  ,  pg 

y  =  -  2u)(x  sin  l)  -  u)  sin  ^[(r  +  z)  cos  1  -  y  sin  t]  -  —  Vy  -  g  — 

o  2  g  r^ 

2  i  PS_  z+r 

z  =  2(u(x  cos  +  u>  cos  -t[(r  +  z)  cos  L  -  y  sin  ^  — —  Vz  -  g  - , 

°  2  P  ^ 

where  x,  y,  z  are  the  radar-centered  cartesian  coordinates  (ft)  of  the 
ballistic  missile,  and 


(1) 


p  =  atmospheric  density  (slug/ft  ) 
g  =  ballistic  coefficient  (lb/ft2) 

v  =  V  x  +  y  +  z 

=  missile  velocity  (ft/sec) 
r2 

g  =  g  —2. 

B  o  2 

r 

1  2 
=  gravitational  acceleration  acting  upon  the  missile  (ft/sec  ) 

2 

g^  =  gravitational  acceleration  at  sea  level  (ft/sec  ) 

r  =  earth  radius  (ft) 
o 

./  2  ”  2  7  ,  i  rF 

r^  =  Y  x  +  y  +  (z  +  rQ) 

=  distance  from  center  of  the  earth  to  the  missile  (ft) 
u)  =  rotation  rate  of  the  earth  (rad/sec) 
t  =  latitude  of  the  radar  site. 
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A  derivation  of  these  equations  is  presented  in  Memorandum  37  and,  hence, 
will  not  be  repeated  here. 


3 


The  x  and  y  axes  point  east  and  north,  respectively,  and  define  a  plane 
tangent  to  the  earth  at  the  radar  site;  the  z  axis  is  perpendicular  to 

this  plane  and  gives  the  vertical  direction  at  the  radar  site.  The 

density  p  is  a  function  of  environmental  factors  and  the  altitude  h, 

which  is  given  by 

h  =  r  -  r *  [x2  +  y2  +  (z  +  r  )2]^  -  r  (2) 

1  o  o  o 


The  ballistic  coefficient  g  is  a  function  of  altitude,  Mach  number, 
and  unknown  reentry  vehicle  parameters. 

Finally,  it  should  be  noted  that  the  above  model  assumes  a  rotating, 

spherical  earth  with  the  ballistic  missile  taken  to  be  a  point  mass. 

That  is,  the  model  given  by  Eqs,  (1) includes  coriolis  terms  (those  con- 

2 

taining  u.0,  centripetal  terms  (those  containing  <n  ),  gravity  gradient 
terms  (those  containing  g),  and  drag  terms  (those  containing  g) . 

Because  g  is  highly  dependent  on  the  characteristics  of  the 
missile  being  tracked,  it  must  be  estimated  together  with  the  position 
and  velocity  state  variables.  This  estimate  is  also  necessary  for  pre¬ 
diction  of  the  future  missile  trajectory.  As  shown  in  Ref.  1,  the  per¬ 
formance  of  the  estimator  is  significantly  improved  if  the  quantity 
p/g  is  estimated,  where  g  is  assumed  to  be  locally  constant  and  p 
is  assumed  to  be  a  locally  exponential  function  of  the  altitude  h,  i.e., 
it  is  assumed  that  p  has  an  altitude  gradient, 


M  _ 


-  K  p  , 


(3) 


in  which  K  is  a  specified  constant  over  each  of  several  ranges  of 
altitude  (h).  Hence,  the  differential  equation  that  models  the  behavior 
of  the  additional  state  variable  p/R  is  given  by 


x  x  +  y  y  +  (z 


=  -  K 


+  r  )z 
o 


(4) 


where  h  is  obtained  from  Eq.  (2) . 
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The  state  vector  x,  which  describes  the  trajectory  of  the  ballistic 
missile,  thus  has  seven  components: 


x3  =  z  xg  =  z 


The  nonlinear  system  equations  (1)  and  (4)  can  then  be  written  in  state- 
variable  form  as  follows: 


X  =  f(x)  = 


2  A 

-2u>(x  cos  l  -  x  sin  ()  +  in  k.  -  i  x  g  Vx  -  g  — 
b  5  x  /  o  **  r 

2  X2 
-2uj(x4  sin  l)  -  (u  sin  li  (ro+x3>cos  l  -  xg  sin  l]  -  ^x^VXg  -  g  — 

x  +r* 

2  ^3  *  q 

2u)(x  cos  l)+  w  cos  (r  +x„)cos  l  -  xn  sin  l]  -  $x  g  Vx  -  g  - 

4  o3  2  7  o  6  r. 


-  K  x„ 


x,  x  +  x  x_  +  (x„  +  r  )x„ 
1  4  2  5  3  o  6 


where  f(x)  is  a  seven-dimensional  vector  function  and 


V=*/ 


s - 5 - T 

X.  +  X+  X. 

4  5  6 


ri  -Vx?  +  x2  +  (X3  +  ro)2' 

The  nonlinear  differential  equation  (6)  is  inexact  because  p  and 
p  do  not  exactly  satisfy  the  assumptions  made  above,  and  also  because 
the  reentry  vehicle  has  additional  degrees  of  freedom  that  have  been 
neglected.  The  actual  trajectory  of  the  ballistic  missile  (i.e.,  the 
trajectory  to  be  estimated)  is  generated  by  a  set  of  equations  that  are 
more  accurate  than  the  model  described  in  Eq.  (6).  In  order  to  account 
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for  these  inaccuracies  as  well  as  actual  system  disturbances,  a  random 
forcing  function  is  introduced  into  Eq .  (5),  yielding 

x  =  f(x>  +  w'  ,  (?) 

where  the  seven-dimensional  vector  w'  is  assumed  to  be  a  zero-mean, 
white  Gaussian  noise  process. 

B.  Radar  Measurement  Equations 

Measurements  of  the  target  vehicle’s  position  are  made  every 
seconds  by  means  of  a  phased-array  radar.  These  measurements  are  made 
of  the  angles  ct>  Y>  and  the  range  r,  which  are  defined  in  the  radar- 
face  coordinate  system  of  Fig.  1.  The  radar-face  coordinates  x  ',  x ',  x' 
are  related  to  the  ground  coordinates  ,  x^,  (or  x,  y,  z)  by  the 
following  two  rotations:  a  rotation  of  9CP  -  t  from  the  vertical  (x3> , 
and  a  rotation  of  £  from  north  (x^) .  Thus, 


*3 


FIG.  1  RADAR  MEASUREMENTS 
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(8) 


“V 

X1 

=  u 

xo 

*3 

u 

_X3 

where 


M  * 


cos  s 
-sin  ? 


0 


sin  -  sin  c 
sin  t  cos  f 
-cos  r 


cos  t  sin  f 
cos  r  cos  " 
sin  t 


and  from  Fig.  1, 


x^  ~  r  sin  a 
x'  a  r  sin  y 

x'  =  r  ll  -  sin“  ;*  -  sin“  y]  ^ 

*3 


(9) 


(10) 


Hence,  the  position  coordinates  x  ,  x  ,  x,  are  related  to  the 
radar  coordinates  r,  a,  y  by  Eqs.  (8),  (9),  and  (10),  with  t  and  P 
being  known  constants  determined  by  the  radar  configuration.  These  non¬ 
linear  relations  can  be  expressed  concisely  as 


r 

a 

Y 


-(xr  *2’  x35  “ 


(11) 


where  l»(x)  is  a  three-dimensional  vector  function.  Thus,  the  nonlinear 
measurement  equation  is  given  by 


z=h(x)+v,  (12) 

where  z  is  a  three-dimensional  vector  containing  the  noise-corrupted 
measurements  of  r,  a,  y  and  the  radar  measurement  noise  v  is  assumed 
to  be  a  zero-mean,  white  Gaussian  noise  process. 
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Ill  EXTENDED  KALMAN  FILTER 


A.  Estimation  Equations 

The  Kalman  filter  equations  specify  an  estimate  of  the  3tate  of  a 
linear,  time-varying  dynamical  system  observed  sequentially  in  the 
presence  of  additive  white  Gaussian  noise;  the  state  estimate  obtained 
at  each  time  is  the  maximum  likelihood  estimate  conditioned  on  all 
measurements  made  up  to  that  time.4  The  Kalman  theory  cannot  be  applied 
directly  to  estimating  the  state  of  a  ballistic  reentry  vehicle  on  the 
basis  of  noise-corrupted  radar  measurements,  because  the  equations  of 
motion  (7)  governing  the  missile  and  the  measurement  equation  (12)  are 
nonlinear  functions  of  the  state  x.  However,  by  first  making  the  ap¬ 
propriate  linearizations,  the  Kalman  theory  can  be  applied  to  the  non¬ 
linear  case.  The  linearization  equations  together  with  the  basic  Kalman 
filter  equations  constitute  what  is  referred  to  as  the  Extended  Kalman 
Filter.  As  shown  in  Refs.  1  and  2,  this  approach  has  been  successfully 
applied  to  trajectory  estimation  problems  in  the  NIKE-X  system.  Justi¬ 
fication  for  this  approach  and  a  derivation  of  the  extended  Kalman  filter 
is  presented  in  Ref.  5. 

Since  the  extended  Kalman  filter  is  to  be  implemented  on  a  digital 
computer,  it  will  be  described  in  its  discrete-time  form  in  this 
memorandum.  The  discrete  time  k  will  refer  to  the  actual  time  k  At, 
where  At  is  the  sampling  interval  (data  processing  interval).  Through¬ 
out  this  report,  x(i/j)  will  denote  the  estimate  of  the  state  x  at 
time  i,  conditioned  on  the  measurements  through  time  j,  and  P(i/j) 
will  denote  the  covariance  of  the  error  in  this  estimate;  i.e,, 

x(i/j)  =  E[x(i)/z(j) ,z(j-i) ,...,z(l) ;x(0/0) ,  P(0/0) ] 

P(i/j)  =  E[ f x(i) -x(i/j) If x(i)-x(i/j) )T/z(j) ,z(j-l) , . . . ,z(l) ;x(0/0) ,P(0/0) j , 

(13) 

where  x(0/0)  is  the  a  priori  state  estimate  and  P(0/0)  is  its  error 
covariance. 

The  measurement  noise  v  given  in  Eq.  (12)  is  a  white  Gaussian 


random  variable  with 


E[v(k)  v(k)  ]  =  R(k)  .  (14) 

The  continuous  random  disturbance  w*  given  in  Eq.  (7)  will  be  replaced 
by  an  equivalent  discrete  random  disturbance  w,  which  is  a  white 
Gaussian  noise  process  with 

E[w(k)]  =  0 

E[w(k)  wT(k)]  =  Q(k)  .  (15a) 

The  value  of  Q(k)  is  a  function  not  only  of  the  time  k,  but  also  of 
the  sampling  interval  At.  If  the  continuous  noise  w'  has  a  covariance 

fn 

given  by  E[w'(t)  w'  (a)]  =  Q'(t)  5(t  -  a),  where  5(0  is  the  Dirac 

delta  function,  then  in  the  discrete-time  form  of  the  missile's  dynamical 

model,  the  discrete  noise  w  will  have  (if  this  model  is  linear)  a  co- 
variance  given  by 

Q(k)  =Q/(k)  At  .  (15b) 

The  initial  state  x(0)  is  assumed  to  be  a  Gaussian  random  variable 
with 

E[x(0) ]  =  5(0/0) 

E[{x(0)  -  x(0/0) }{x(0)  -  x(0/0) }T]  =  P(0/0)  .  (16) 

Furthermore,  it  is  assumed  that  v(k) ,  w(k) ,  and  x(0)  are  mutually 
uncorrelated. 

The  resulting  estimation  equations,  which  are  given  by  the  extended 
Kalman  filter,  can  be  considered  as  consisting  of  two  parts:  prediction 
and  correction. 

1.  Predict  ion 

Given  x(k/k) ,  the  estimate  of  the  state  at  time  k,  and  the 

assumption  that  E[w(k)]  =  0,  the  predicted  state  can  be  obtained  from 

* 

the  nonlinear  equations  of  motion  (7): 

x(k+l/k)  =  x(k/k)  +  f[x(k/k)]At  ,  (17) 

with  the  covariance  of  the  error  in  the  state  prediction  x(k+l/k)  given  by 

More  accurate  integration  formulas  can  be  used  if  necessary. 
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(18) 


P(k+l/k)  =  §  (k+1  ,k)  P(k/k)  $T(kn,k)  +  Q'(k)  At, 

where  §(k+l,k)  is  obtained  by  linearization  of  Eq.  (7)  about  the 
state  estimate  x(k/k) — i.e., 


§(k+l,k)  =  I  +  F[x(k/k)]At  , 
in  which  F  is  the  7x7  matrix  with  components 

Sf . (x) 


V*  - 


(19) 


(20) 


From  Eq.  (6)  , 
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where 


2  -1  ,  „  2-3 

F41  -  8  rl  *  3  S  \ri 


F  „  =  3  g  x,  x  r 
42  b  1  2  1 
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,-l 
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F67  =  -  4  V  *6  g0 
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P76  *  '  *7  K  T'l  (*3  +  'o' 

=  -  K  r^1  lx,  x4  +  x2  +  (*3  *  ro)x6!  . 


2.  Correction 

The  state  prediction  x(k+l/k)  is  then  corrected  by  using 
z(k+l)  ,  the  actual  measurement  at  time  k+1,  and  z(k+l/k) ,  the  pre¬ 
dicted  measurement  for  time  k+1,  which  is  obtained  from  the  nonlinear 


measurement  equation  (12)  by  using  the  assumption  that  E[v(k+1)]  =  0: 

z(k+l/k)  =  h[x(k+l/k) 3  «  (22) 

Hence,  the  estimate  of  the  state  at  time  k+1  is  given  by 
x(k+l/k+l)=  x(k+l/k)  +  W(k+l)[z(k+l)  -  2(k+l/k) ]  ,  (23) 

where  the  weighting  matrix  is 

W(k+1)  =  P(k+l/k)HT(k+l)fR(k+l)  +  H(k+l)P(k+l/k)HT(k+l) ]_1  (24) 

and  H(k+1)  is  obtained  by  linearization  of  Eq.  (12)  about  the  state 
prediction  x(k+l/k) — i.e., 

H(k+1)  =  H[x(k+l/k) ]  ,  (25) 


in  which  H  is  the  3x"  matrix  with  components 


dh  (x) 

HU(i>  ■  1^- 


(26) 


From  Eqs.  (8)  through  (11), 


ah." 

Hn 
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0 

0 

0 

0 

H  = 

l 

= 

H21 

H22 

H23 

0 

0 

0 

0 

H31 

H32 

H33 

0 

0 

0 

0 

(27) 


where  the  H  are  most  conveniently  written  in  terms  of  the  radar  face 
coordinates  x',  x^ ,  x^: 
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X  x  * 

Hn  =  ~  cos  §  +  ~  Sin  7  sin  § 


+  —  cos  t  sin  5 
r 


M12  *  [‘’'s'2  +  <’t3>2  ]  «»  I 
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x'  , 

”  L  2)2  +  ^xp2  J  cos  t  sin  £ 
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31  r  r 

x'  . 
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=  ~  [|X1^2  +  ^X3^2J  C°S  T 


"  IT'  *  <X3>2  J!  sin  T 
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For  use  in  Eq,  (27),  the  values  of  x' ,  x',  x'  are  determined  from  the 

*  W  W  ^ 

state  prediction  x(k+l/k)  by  the  linear  transformation  M  ,  where  M 
is  given  in  Eq.  (9): 

x'l  ^^(k/k-l) 

*2  =  M**1  S2(k/k-1> 

x'  *3(k/k-l) 

The  covariance  of  the  error  in  the  estimate  x(k+l/k+l)  is 

POd-l/k+1)  =  (I  -  W(k+l)H(k+l)  ]P(k+l/k) 

=  P(k+l/k)~P(k+l/k)HT(k+l)[R(k+l)+H(k+l)P(k+l/k)HT(k+l)]-1 

•  H(k+l)P(k+l/k)  .  (29) 

Equations  (17),  (18),  (23),  (24),  and  (29),  together  with  the 
matrices  defined  in  Eqs.  (19)  and  (25)  ,  comprise  the  extended  Kalman 
filter.  The  a  priori  state  estimate  x(0/0)  and  its  error  covariance 
P(0/0)  are  used  to  initialize  these  recursive  equations.  The  measure¬ 
ment  noise  covariance  R  is  specified  by  the  characteristics  of  the 
phased-array  radar  model.  The  random  disturbance  covariance  Q  is 
determined  in  order  to  compensate  for  inaccuracies  in  the  model  of  the 
equations  of  motion  (see  Secs.  V  and  VI). 

It  should  be  noted  that  in  the  extended  Kalman  filter,  the 
nonlinear  equations  (7)  and  (12)  are  used  to  obtain  the  state  prediction 
(Eq.  (17)]  and  the  predicted  measurement  [Eq.  (22)],  respectively.  The 
linearization  of  Eqs,  (7)  and  (12)  ,  used  to  obtain  §  and  H,  is  only 
employed  to  calculate  the  covariances  and  the  weighting  matrix  in 
Eqs.  (18),  (24),  and  (29). 

B.  Basic  Flow  Diagram 

The  basic  equations  of  the  extended  Kalman  filter  (see  Sec.  Ill  A) 

can  be  displayed  in  the  form  of  a  flow  diagram,  which  is  a  simplified 

pictorial  representation  of  the  digital  computer  program  that  has  been 

developed.  This  is  shown  in  Fig.  2,  where  k  is  the  maximum  number 

max 

of  times  that  the  filter  processes  the  measurements.  Thus,  there  are 
seven  major  computations  involved  in  each  filtering  iteration — as  in¬ 
dicated  by  the  seven  numbered  blocks  in  Fig.  2  (these  will  be  referred 
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FIG.  2  FLOW  DIAGRAM  FOR  EXTENDED  KALMAN  FILTER 
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to  as  Steps  1  through  7) .  Since  the  state  vector  of  the  ballistic 
missile  has  seven  components,  some  of  the  calculations  shown  in  Fig.  2 
require  considerable  computation  time. 

In  the  next  section,  the  contribution  made  by  each  of  the  seven 
basic  calculations  to  the  total  computational  time  for  one  filtering 
Iteration  will  be  found.  These  timing  estimates  provide  the  basis  for 
making  modifications  in  the  filter  design  to  decrease  the  average  compu¬ 
tation  time  per  iteration. 

C.  Computation  Time  Considerations 

In  a  previous  study,1  when  the  extended  Kalman  filter  equations 
were  implemented  on  the  IBM  7090  for  endoatmcspheric  trajectory  estima¬ 
tion,  the  filter  ran  approximately  6  times  slower  than  real  time  with  a 
data  rate  of  20  measurements  per  second  (for  a  single  target);  i.e.,  it 
took  6  seconds  of  computer  time  to  process  the  measurements  obtained  in 
1  second.  Since  the  speeo  of  the  Kalman  filter  is  computer-dependent, 
one  can  improve  the  speed  with  a  faster  computer.  In  fact,  with  the  same 
data  rate  of  20  measurements  per  second,  the  Univac  1108  runs  approxi¬ 
mately  5  times  faster  than  real  time  (for  a  single  target)  for  endo- 
atmospheric  estimation.  Even  at  five  times  faster  than  real  time  (for 
a  single  target) ,  the  speed  of  the  extended  Kalman  filter  may  still  be 
inadequate  for  some  situations.  For  instance,  if  n  missiles  must  be 
tracked  and  if  the  speed  of  the  Kalman  filter  is  5  times  real  time, 
then  for  n>  5,  the  filter  cannot  compute  state  estimates  for  all  of 
the  n  missiles  in  real  time  (at  a  data  rate  of  20  measurements  per 
second) . 

For  exoatmospheric  trajectory  estimation  (for  which  the  filter  has 
six  states  instead  of  seven  and  is  programmed  in  double  precision — see 
Sec.  VI) ,  the  Univac  1108  runs  approximately  100  times  faster  than  real 
time  (for  a  single  target)  with  a  data  rate  of  one  measurement  per 
second . 

Thus,  particularly  in  a  multi-threat  situation  for  endoatmospheric 
estimation,  it  is  essential  to  improve  the  computational  efficiency  of 
the  extended  Kalman  filter.  This  should  be  done  with  a  minimal  amount 
of  degradation  in  filter  performance.  For  ballistic-trajectory 
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estimation,  where  the  dynamical  description  is  given  by  Eqs.  (6),  (7), 
and  (12) ,  we  know  the  particular  structure  of  the  linearized  differential 
equation  matrix  F,  the  transition  matrix  §,  the  linearized  measurement 
matrix  H,  the  random  disturbance  covariance  Q,  and  the  measurement- 
noise  covariance  R.  Therefore,  we  can  take  advantage  of  this  knowledge 
by  writing  a  computer  program,  which  reproduces  the  flow  diagram  (Fig.  2) 
in  a  computationally  efficient  manner. 

In  the  extended  Kalman  filter  program,  which  has  been  developed, 

the  seven  major  computations  per  iteration  given  in  Fig.  2  involve 

*• 

various  numerical  operations  as  listed  in  Table  I. 


Table! 

ARITHMETIC  OPERATIONS  FOR  THE  EXTENDED  KALMAN  FILTER 
(ENDOATMOSPHERIC  CASES) 


COMPUTATIONAL 

OPERATION 

(see  Fig.  2) 

NUMBER  OF 
MULTIPLICATIONS 
OR 

DIVISIONS 

CALCULATE  LINEARIZED  DIFFERENTIAL- 
EQUATION  MATRIX,  EQ.  (21) 

89 

CALCULATE  TRANSITION  MATRIX, 

EQ.  (19) 

28 

CALCULATE  PREDICTED  STATE, 

EQ.(17) 

30 

CALCULATE  ERROR  COVARIANCE 

IN  PREDICTED  STATE,  EQ.  (18) 

693 

CALCULATE  WEIGHTING  MATRIX, 

EQS.  (24)  AND  (27) 

102 

CALCULATE  STATE  ESTIMATE, 

EQ.  (23) 

21 

CALCULATE  ERROR  COVARIANCE  IN 
STATE  ESTIMATE,  EQ.  (29) 

147 

NUMBER  OF 
ADDITIONS 
OR 

SUBTRACTIONS 


NUMBER  OF 
SQUARE  ROOTS 


TOTALS 


1110 


100? 


2 


The  enumeration  of  arithmetic  operations  in  Table  I  gives  an  indication 
of  the  distribution  of  the  total  computation  time  per  iteration  among 
the  seven  basic  calculations.  For  instance,  if  it  can  be  assumed  that 
the  operation  of  multiplication  (or  division)  requires  a  computation 
time  of  an  order  of  magnitude  larger  than  that  for  an  addition  (or 
subtraction) ,  then  the  calculation  of  the  predicted  error  covariance 
(Step  4)  takes  about  63%  of  the  total  computation  time.  Clearly,  the 
error  covariance  matrix  calculations  of  Steps  4  and  7  are  the  most 
involved. 


The  next  section  discusses  various  approaches  that  may  be  used  to 
modify  the  Kalman  filtering  algorithm  in  order  to  reduce  the  computational 
requirements.  Particular  emphasis  is  placed  upon  simplifying  and  modi¬ 
fying  the  two  error  covariance  calculations,  since  they  are  the  most 
time-consuming  calculations  in  the  Kalman  filter. 


Approaches  to  Reducing  the  Computational  Requirements 


One  approach  to  reducing  the  computational  requirements  of  the 
Kalman  filter  was  discussed  in  detail  in  Ref.  1.  There  an  approximation 
to  the  weighting  matrix  W(k)  was  precomputed  off-line,  rather  than 
calculated  recursively  on-line.  The  weighting  matrix  values  for  repre¬ 
sentative  trajectories  are  then  stored  in  the  computer  memory.  The 
precomputed  weighting  matrix  technique  has  some  disadvantages.  First, 
the  weighting  matrix  W(k)  not  only  is  a  function  of  the  iteration  time 
k,  but  it  also  is  altitude-dependent.  Second,  W(k) ,  even  for  the  same 
altitude,  is  dependent  to  some  degree  upon  the  geometry  of  the  trajectory 
and  the  ballistic  coefficient  (at  lower  altitudes).  Consequently,  for 
accurate  estimation,  many  representative  weighting-matrix  histories 
should  be  stored  in  memory  and  employed  accordingly.  The  computational 
advantages  of  the  Kalman  filter  using  a  precomputed  weighting  matrix  are 
impressive,  however.  In  fact,  for  endoatmospheric  estimation,  the  speed 
of  this  filter  is  approximately  15  times  that  of  the  fully  implemented 
Kalman  filter. 


Another  approach  for  improving  the  real-time  capability  of  the 
Kalman  filter  is  to  decrease  the  data  rate  until  it  falls  within  the 
filter's  real-time  capability.  One  disadvantage  of  this  approach  is  that 


the  estimation  is  based  upon  less  data  and,  thus,  is  inherently  less 
accurate.  As  At,  the  time  between  measurements,  is  increased,  the 
calculations  yielding  $(k+l,k)  and  x(k+l/k)  should  be  refined  so 
that  their  accuracy  is  not  impaired;  these  refinements  consume  some  of 
the  computational  speed  gained  by  decreasing  the  data  rate.  Additionally, 
the  linearization  assumptions  for  the  extended  Kalman  filter  may  no 
longer  be  valid  for  large  values  of  At .  Clearly  the  speed  of  the  Kal- 
men  filter  relative  to  real  time  is  proportional  to  At,  so  that  if  the 
filter  is  6  times  slower  than  real  time  for  a  data  rate  of  20  measure¬ 
ments  per  second  (At  =  0.05  sec),  then  it  will  be  faster  than  real  time 
only  for  data  rates  less  than  3.33  measurements  per  second.  In  Ref.  1, 
it  was  shown  that  the  performance  of  the  extended  Kalman  filter  for 
endoatmospheric  estimation  was  not  markedly  affected  when  the  data  rate 
was  decreased  by  up  to  a  factor  of  5  from  a  basic  data  rate  of  20 
measurements  per  second.  Some  further  results  on  this  approach  will  be 
presented  in  Secs.  V  and  VI. 

IV  PIECEV/ISE-RECURSIVE  KALMAN  FILTER 

In  this  study  we  have  combined  the  advantages  of  both  of  the 

approaches  discussed  in  See.  Ill  D.  The  filter  algorithm  considered 

here  calculates  the  weighting  matrix  W(k+1)  every  AT  seconds,  and 

uses  this  weighting  matrix  for  calculating  the  state  estimate  x(i/i) 

At 

at  intervals  of  At  sec  for  i  =  k+l,  k+2 . k+N,  where  —  =  N  £  1. 

This  is  equivalent  to  calculating  W(k+1)  at  time  k+l  and  retaining 
this  value  of  the  weighting  matrix  (to  be  used  in  computing  the  state 
estimate)  for  the  times  k+l,  k+2,  ....  k+N,  and  then  calculating 
w[ k+(N+l) ]  anew.  The  design  parameters  of  the  filter  are  then  At,  the 
sample  interval  (data  processing  interval),  and  AT ,  the  interval  between 
weighting-matrix  calculations.  This  filter  has  a  "piecewise-precomputed" 
weighting  matrix  and  its  data  rate  may  also  be  varied.  Thus,  the  computa¬ 
tional  advantages  of  both  precomputation  and  data-rate  reduction  can  be 
achieved.  This  algorithm  also  avoids  the  off-line  nature  of  a  Kalman  filter 
that  uses  a  precomputed  approximation  to  the  weighting  matrix;  in  this 
approach,  the  weighting-matrix  calculations  are  performed  recursively  at 
Intervals  of  At  seconds.  This  computational  algorithm,  which  will  be 
referred  to  as  the  Plecewise-Recurslve  (Extended)  Kalman  Filter,  will  be 
described  more  completely  below. 
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The  modified  Kalman  filter  algorithm  described  above  operates  as 
follows:  In  the  first  cycle,  the  equations  for  the  state  prediction 
and  estimate  are  computed  just  as  they  are  for  the  fully  implemented 
Kalman  filter  [see  Sec.  Ill  A,  Eqs.  (17),  (22),  and  (23)], 

x(k+l/k)  =  x(kA)  +  flx(k/k)Mt  (30) 

x(k+l/k+l)  =  x(k+l/k)  +  W(l)[z(k+1)  -  z(k+l/k) j  (31) 

k  *  0,  1 . N-i 

The  weighting  matrix  W(l)  for  these  K  iterations  is  obtained  from 
Eqs.  (18)  and  (24), 

V(l)  =  P(1/0)HT(1)[R(D  +  H(l)  P(1/0)HT(1)]"L  ,  (32) 

where 

P(l/0)  =  §(1,0)  P(0/0)  $T(l,0)  +  Q'(0Ut  (33) 

and 

§(1,0)  =  I  +  F[x(0/0)]At 

The  calculations  in  Steps  1,  2,  3,  4,  5,  and  7  of  Fig.  2  are  not  done 
for  k  =  1,  2,  ... ,  N-l. 

Before  proceeding  to  the  next  cycle,  this  algorithm  requires  that 

P(N/N) ,  the  covariance  of  the  error  in  the  state  estimate  x(N/N)  at 

the  end  of  the  previous  cycle,  be  computed.  As  shown  in  Ref.  6,  this 

* 

covariance  can  be  obtained  by  solving  the  recursive  relations: 

P(k-tl/k^l)  =  [I  -  W(l)H(kxl)  ]  P(kn/k)  [i  -  W(l)H(k+l)  ] 

x  W ( 1)  R ( k+1) \VT ( 1)  (34) 

P(k+l/k)  =  i (k+1 ,k)P(k/k)§T(k+l ,k)  +Q'(k)At  (35) 

k  =  0,  1 . N-l 

However,  this  calculation  of  P(N/N)  defeats  the  very  purpose  of  the 
piecewise-recursi ve  Kalman  filter — namely,  to  reduce  the  computation 


These  relations  exactly  define  the  propagation  of  the  covariance  for 
either  optimal  or  sub-optimal  values  of  the  weighting  matrix. 
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time.  Therefore,  simplifications  for  or  approximations  to  Eqs.  (34)  and 
(35)  are  essential. 


The  approximation  made  here  is  that  we  can  consider  the  filtering 
done  at  the  times  k  =  1,  2,  ...,  N  as  constituting  a  presmoothing  of 
measured  data.  In  effect,  it  can  be  assumed  that  at  the  end  of  the  time 
interval  [0,  At],  where  k  =  N,  one  measurement  (with  an  appropriately 
reduced  noise  covariance)  is  processed.  It  can  be  shown  that  for  the 
most  optimistic  conditions  this  reduced  measurement  noise  covariance  is 
given  by 


R*  (N) 


R(N) 

N 


(36) 


where  R(N)  is  the  average  value  of  the  measurement -noise  covariances 

for  the  times  k  =  1,  2,  ....  N  (At  =  NAt) ;  for  practical  purposes, 

_  * 

R(N)  is  taken  to  be  equal  to  R(N).  It  is  asserted  that  R  (N) 

represents  (approximately)  the  effective  measurement  noise  over  the 

interval  [o,  At],  if  all  N  measurements  are  processed  by  the  Kalman 

filter. 


Using  this  reduced  measurement-noise  covariance,  P(N/N)  is  then 
determined  [see  Sec.  Ill  A,  Eqs.  (18)  and  (29)]  from  the  equations 

P(N/N)  =  P(N/0)  -  P(N/0)  HT(N)  +  H(N)  P(N/0)  HT(N)]_1 

•  H(N)  P(N/0)  ,  07) 

where 

P(N/0)  =  $(N,0)  P(0/0)  $T(N,0)  +  Q'(0)  NAt  (38) 


and 

5(N,0)  =  I  +  F[x(0/0)]NAt 


This  approach  avoids  the  recursive  solution  of  P(N/W  and  has  been 
found  to  yield  satisfactory  filter  performance  as  shown  in  Secs.  V  and 
VI. 

At  k  =  N,  the  cycle  is  repeated;  i.e.,  the  state  prediction  and 
estimate  are  given  by 

x(k+l/k)  =  x(k/k)  +f[x(k/k)]At  (39) 
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x(k+l/k+l)  =  x(k+l/k)  +  W(N+l)[z(k+l)  -  z(k+l/k) ]  (40) 

k  =  N,  N+l . 2N-1 

The  weighting  matrix  W(N+1)  for  these  N  iterations  is  obtained  as 

W(N+1)  =  P(N+i/N)  HT(N+1) [r(N+1)  +  H(N+1)  P(N+1/N)  HT(N+1)]_1,  (41> 

where 

P(N+1/N)  =  § (N+l ,N)  P(N/N)  ST(N+1,N)  +  Q/(N)At  (42) 

and 

§ (N+l ,N)  =  I  +  F(x(N/N) ] At 

Again,  the  computations  in  Steps  1,  2,  4,  5,  and  7  in  Fig.  2  are  omitted 
for  k  =  N+l,  N+2 ,  ...,  2N-1.  Next,  P(2N/2N) ,  the  covariance  of  the 
error  in  the  state  estimate  x(2N/2N)  at  the  end  of  this  cycle,  is  cal¬ 
culated  from 

P(2N/2N)  =  P(2N/N)  -  P(2N/N)HT(2N)  +  H(2N)P(2N/N)HT(2N)  ]_1 

•  H(2N)  P(2N/N)  ,  (43) 

where 

P(2N/N)  =  $(2N,N)  P(N/N)$T(2N,N)  +  Q '(N)NAt  (44) 

and 

$ (2N,N)  =  I  +  F[x(N/N)]NAt 

Whenever  k  is  a  multiple  of  N,  this  cycle  is  repeated  in  the 
manner  outlined  above. 

The  use  of  the  piecewise-recursive  Kalman  filter  results  in  a 
reduction  of  computational  requirements  whenever  N  >  2  (i.e., 

At  >  2  At).  The  reason  for  this  is  that  the  modified  algorithm  nearly 
doubles  the  number  of  numerical  calculations,  with  respect  to  the 
basic  Kalman  filter  presented  in  Fig.  2,  whenever  k  is  a  multiple  of 
N;  but  for  the  next  N-l  sample  points,  the  filter's  weighting  matrix 
is  fixed  so  that  the  number  of  numerical  calculations  is  very  small  for 
these  iterations.  Thus,  there  is  no  computational  advantage  of  this 
algorithm  unless  N  >  2. 
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The  expected  decrease  in  the  computational  time  for  the  piecewise- 

recursive  Kalman  filter  is  a  function  of  the  size  of  At  relative  to 

At,  as  well  as  the  value  of  At.  Let  S  be  the  ratio  of  the  compu- 

o 

tational  speed  of  a  Kalman  filter  using  a  precomputed  weighting  matrix 
to  the  computational  speed  of  the  fully  implemented  Kalman  filter  at 
the  same  data  rate  (Sq  is  independent  of  the  sample  interval  At) . 

Since  the  piecewise-recursive  Kalman  filter  described  above  is  equivalent 
to  using  a  precomputed  weighting  matrix  for  N-l  iterations  (where 
N  =  At /At)  before  a  new  weighting  matrix  is  computed,  this  filter  is  S 

o 

times  faster  than  the  fully  implemented  Kalman  filter  for  these  N-l 

iterations.  As  shown  above,  for  the  first  iteration  of  each  cycle  of 

N  iterations,  the  piecewise-recursive  Kalman  filter  is  approximately 

twice  as  slow  as  the  fully  implemented  Kalman  filter.  Thus ,  if  At  is 

K 

the  computational  time  consumed  by  the  fully  implemented  Kalman  filter 

for  one  iteration,  then  the  time  consumed  by  the  piecewise-recursive 

Kalman  filter  is  approximately  2At  if  the  iteration  is  the  first  in 

the  cycle  of  N  iterations,  and  At„/S  if  the  iteration  is  one  of 

K  O 

the  remaining  N-l  iterations  before  a  new  weighting  matrix  is  calcu¬ 
lated. 


Since  the  total  computational  time  for  N  data  points  with  the 

fully  implemented  Kalman  filter  is  NAt  ,  and  with  the  piecewise-recursive 

K 

Kalman  filter  it  is  2At„  +  (N-l)  At„/S  ,  it  follows  that  the  ratio  of 

K  K  o 

their  average  iteration  times  is  approximately 


2At„  +  (N-l)  At„/S 

K _ K  O 

N  AtK 


(45) 


Thus,  the  speed  of  the  piecewise-recursive  Kalman  filter,  relative  to 
the  fully  implemented  Kalman  filter,  is  approximately 


N  S 

s  -  .  -2 _ 

N  "  2S  +  (N-l) 
o 


(46) 


If  the  speed  of  the  fully  implemented  Kalman  filter  relative  to 
real  time  (when  the  sample  interval  is  At  )  is  S  ,  then  its  speed 

O  K 

relative  to  real  time  when  the  sample  interval  is  At  will  be 
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(47) 


Consequently,  the  speed  of  the  piecewise-recursive  Kalmar,  filter  with 
respect  to  real  time  will  be  approximately 


R  vAt 


2S  +  (N-l) 
o 


As  an  illustration  of  the  expected  computational  advantages  in 

endoatmospheric  estimation  of  such  a  filter  design  using  the  Univac  1108 

at  a  sampling  interval  of  At  =  0.05  sec,  the  fully  implemented  Kalman 

o 

filter  is  about  5  times  faster  than  real  time  (i.e.,  S  =  5)  ;  Table  I 

R 

indicates  that  the  speed  of  a  Kalman  filter  using  a  precomputed  weighting 

matrix  relative  to  the  fully  implemented  Kalman  filter  is  about  15 

(i.e.,  S  =  15).  Incidentally,  it  should  be  noted  that  S  is  computer- 
O  R 

dependent,  whereas  S  is  not.  Thus,  for  the  endoatmospheric  estimation 

o 


on  the  Univac  1108,  3 


is  given  by 


©15N  j 

30  +  (N-l) 

*—  -J 


Note  that  S  ,  the  speed  of  the  piecewise-recursive  Kalman  filter 
fl  ?  fat 

relative  to  real  time,  increases  linearly  with  At;  for  fixed  At  it  is 
bounded  by  the  speed  of  a  Kalman  filter  with  a  precomputed  weighting 
matrix. 

To  obtain  some  measure  of  the  improvement  in  the  computational 

speed  for  the  endoatmospheric  cases,  let  us  calculate,  using  Eq.  (49), 

At 

the  iteration  times  for  various  values  of  N  =  —  .  The  average  time 

iit 

for  one  iteration  on  the  Univac  1108  is  given  by  At/S.r  ;  thus,  for 

N ,  At 

N  =  5,  10 ,  20 ,  and  30,  we  obtain  the  average  iteration  times  of  0.0045, 
0.0026,  0.0016,  and  0.0013  sec,  respectively.  (Recall  that  the  iteration 
time  for  the  fully  implemented  Kalman  filter  is  0.0100  sec.)  These 
estimates  of  the  average  iteration  time  for  the  piecewise-recursive 
Kalman  filter  are  in  very  close  agreement  with  the  actual  values  found 
experimentally  as  will  be  shown  in  Sec.  V.  For  exoatmospheric  estima¬ 
tion,  the  average  iteration  times  of  the  piecewise-recursive  Kalman 

At 

filter  (for  various  values  of  N  =  —  )  are  given  in  Sec.  VI. 

At 
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In  Secs.  V  and  VI,  we  will  verify  the  expected  computational  ad¬ 
vantages  of  the  piecewise-recursive  Kalman  filter  and  investigate  the 
performance  of  this  filter  for  endoatmospheric  and  exoatmospheric  esti¬ 
mation. 

Since  the  piecewise-recursive  Kalman  filter  algorithm  essentially 
employs  a  weighting  matrix  whose  elements  are  piecewise-constant 
functions  of  the  time  k,  one  could  theoretically  determine  the  optimal 
values  of  the  weighting  matrix  elements,  subject  to  this  piecewise- 
constant  constraint.  This  appears  to  be  a  difficult  problem  to  solve 
analytically.  The  piecewise-recursive  Kalman  filter  described  here  can 
be  considered  as  an  approximate  solution  to  this  problem. 

V  N  tUCAL  RESULTS  FOR  ENDOATMOSPHERIC  TRAJECTORY  ESTIMATION 


A.  Description  of  the  Test  Cases  and  Filter  Design  Parameters 

Data  from  two  representative  endoatmospheric  trajectories  generated 
by  simulation  using  an  accurate  model  of  reentry  dynamics6  was  used  to 
evaluate  the  performance  of  the  piecewise-recursive  Kalman  filter  des¬ 
cribed  in  Sec.  IV.  The  computer  program  used  in  generating  these 
trajectories  was  developed  by  M.I.T.  Lincoln  Laboratory. 

The  reentry  angle  for  both  trajectories  is  approximately  22°. 

Case  1  is  a  trajectory  of  a  high-g  missile,  which  is  tracked  starting 
at  an  altitude  of  about  205,000  ft.  This  trajectory  impacts  approximately 
at  the  radar  site.  Case  2  is  a  trajectory  of  a  low-B  missile,  which  is 
tracked  from  an  altitude  of  about  153,000  ft.  This  trajectory  impacts 
about  80,000  ft  east  and  30,000  ft  north  of  the  radar  site.  Cases  1 
and  2  are  essentially  the  same  as  Cases  1  and  4,  respectively,  of  Ref.  1. 

The  filtering  computations  start  with  an  initial  state  estimate 
x(0/0)  and  its  error  covariance  matrix  P(0/0) .  The  initial  state 
estimates  for  the  endoatmospheric  cases,  generated  as  in  Ref.  7,  are 
given  in  Table  II. 

Based  on  the  work  presented  in  Ref.  1,  the  initial  error  covariance 
{for  the  endoatmospheric  trajectories)  is  chosen  to  be  the  7  x  7 
diagonal  matrix 
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P(0/0)  = 


5  X  10 


The  measurement -noise  covariance  R(k)  is  obtained  from  a  model  of 
actual  radar  noise  characteristics  of  the  M5R  based  on  work  by  Bell 
Telephone  Laboratories.8  The  expressions  for  the  elements  of  this  co- 
variance  are  found  in  Ref.  9.  This  radar  model  is  used  to  generate  the 
measurement  noise  of  Eq.  (12)  in  the  simulation.  Observations  may  be 
taken  as  often  as  every  0.05  sec  for  the  endoatmospheric  cases,  but  the 
filter  may  process  fewer  measurements  than  are  taken  (i.e.,  it  &  0.05 

sec)  . 

Toble  II 

INITIAL  CONDITIONS  FOR  ENDOATMOSPHERIC  TRAJECTORIES 


(ft/sec)  (ft.  sec) 


*7 

(lb.  ft ‘) 


338110  338110  199910  15297  15297  8653 

290952  292135  199911  16046  14466  12751 


4.395  x  10-  10 

6.118  x  10- 10 


83  1  4098 


21633  8653  7.950  x  10' 

20891  11592  2.376  x  10" 


41583  68 


742  2939 


The  random  disturbance  covariance  Q(k)  [see  Eqs.  (15a)  and  (15b)], 
compensates  for  the  model  inaccuracies  as  described  in  Sec.  II  A.  Based 
on  studies  described  in  Ref.  1,  this  covariance  (for  the  endoatmospheric 
trajectories)  is  chosen  to  be  the  7  x  7  diagonal  matrix 


(51) 


where 


and 

B. 


3.2X10"19 

,  x3 Si 60000 

-18 

1.0x10 

,95000ex3< 160000 

2.2Xl0“15exp(-5xl0_5x3) 

,25000ex  <  95000 

u 

2.2xl0"lD[l+1.5xl0_3(2.7xl04-x  ) ]exp(-5xl0_5x„) ,14700sx„<  25000 

O  O  o 

3.6xlO~14[l+1.5xl0~J(2.7xlO4-x  ) ]exp(-5xl0“°x  ) ,  x„<  14700 

V.  J  O  o 

xg  is  the  present  vertical  distance  estimate. 

Performance  of  the  Fully  implemented  Extended  Kalman  Filter 


For  the  endoatmospheric  trajectories,  the  standard  which  is  used  to 
evaluate  the  performance  of  the  piecewise-recursive  Kalman  filter  is  the 
fully  implemented  (extended)  Kalman  filter  with  a  sample  interval  of 
At  =  0.05  sec.  The  recursive  equations  for  the  extended  Kalman  filter 
are  given  in  Sec.  Ill  A.  This  filter  was  used  to  estimate  the  ballistic 
trajectories  for  Cases  1  and  2  using  simulated  trajectories  and  measure¬ 
ment  noise.  The  filter  yields  state  estimates  x(k/k)  at  sample 
intervals  of  At  =  0.05  sec. 


The  magnitudes  of  the  measurement  errors  in  position  are  shown  in 
Figs,  3a  and  b  for  Cases  1  and  2.  In  order  to  relate  altitude  and  time, 
these  magnitudes  are  plotted  as  a  function  of  both  variables;  in  all 
other  figures,  the  plots  are  a  function  of  time  only.  Figures  4a,  b, 
and  d  show  the  magnitudes  of  the  estimation  errors  for  position,  velocity, 
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5000 


sec 


and  acceleration,  respectively,  for  Case  1.  Figure  dc  shows  both 
actual  and  estimated  P,  where  the  estimated  B  is  obtained  by  assuming 
§(k/k)  =  o(h(k/k) ]/x7(k/k)  in  which  fi(k/k)  is  the  estimated  altitude 
obtained  from  Eq.  (2).  It  should  be  noted  that  acceleration  is  not 
computed  directly  by  the  filter  but  is  calculated  from  the  state  esti¬ 
mate  using  Eqs.  (1)  of  Sec.  II  A.  Figures  5a,  b,  c,  and  d  show  the 
corresponding  estimation  errors  for  Case  2. 

Note  that  after  3.5  sec  of  filtering,  the  magnitude  of  the  errors 

in  estimated  position  in  both  cases  is  never  any  worse  than  830  ft  (and 

in  Case  1,  the  initial  position  error  is  approximately  65,000  ft). 

After  3.5  sec,  the  magnitude  of  error  in  estimated  velocity  in  both 

cases  is  never  any  worse  than  285  ft /sec.  The  magnitude  of  error  in 

2 

estimated  3  for  Case  1  is  never  more  than  665  Ib/ft  ,  and  for  Case  2 

2 

it  is  never  more  than  105  lb/ft  after  3.5  sec  of  filtering. 

C.  Performance  of  the  Fully  Implemented  Extended  Kalman  Filter  with 

Different  Data  Rates 

The  simplest  scheme  for  improving  the  real-time  capability  of  the 
Kalman  filter  is  simply  to  decrease  the  data  rate.  If  the  sample 
intervals  are  lengthened  and  the  additional  time  between  samples  is  not 
used  to  presmooth  the  measurements  (i.e.,  the  intermediate  measurements 
are  discarded),  then  obviously  the  filter  performance  will  deteriorate 
to  some  extent.  Reference  1  discusses  the  performance  of  the  extended 
Kalman  filter  for  four  representative  endoatmospheric  cases.  However, 
no  attempt  was  made  in  Ref.  1  to  investigate  the  effect  of  both  decreasing 
the  data  rate  and  presmoothing  the  intermediate  measurements.  The  present 
study  considers  the  situation  where  the  presmoothing  is  dene  by  fitting 
the  data  to  a  curve  that  approximates  the  trajectory  of  the  ballistic 
missile.  Of  course,  in  practice,  data  is  presmoothed  by  curve-fitting 
with  polynomials  that  are  not  precise  representations  of  missile  tra¬ 
jectories  . 

The  assumption  of  data  presmoothing  implies  that  in  the  filter 
equations  the  measurement-noise  covariance  should  be  reduced  accordingly. 
If  the  measurement-noise  covariance  at  a  particular  time  is  R° (k)  for 
a  sample  interval  of  AtQ,  then  one  approximation  to  its  reduced  value, 
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FIG  5  ESTIMATION  ERRORS  FOR  FULLY  “IMPLEMENTED  KALMAN  FILTER  —  CASE  2.  \t  0.05  see 


with  a  larger  sample  interval  At  and  with  presmoothing ,  is 
Ato 

R(k)  =  R°(k)  .  To  see  how  this  is  derived,  suppose  that  the  measure¬ 

ment  vector  is  a  linear  function  of  time  t  in  the  interval  [t  , 

1  K 

t  +  At]  with  additive  noise  v^  <i.e.,  =  p  +  st^  +  v^)  ,  and  its 

rate  of  change,  s,  at  t  =  t  is  known  (equivalently,  our  estimate 

K 

of  the  velocity  s  equals  its  actual  value).  Now,  we  estimate  the 

unknown  parameter  p  by  subtracting  st^  from  the  measurements  at  the 

times  t  =  t,  ,  t  +  At  ,  t  +  2At . t  +  At  (where  At  =  n  At  ) 

lkkoko  k  o 

and  average  the  resulting  quantities.  Let  p  designate  the  computed 

average  and  z(t  +  At)  =  p  +  s(t  +  At)  be  the  presmoothed  version  of 
K  K 

the  measurement  z(t,  +  t) ,  where  v.  is  zero  mean, 

k  i 

Noting  that  E[pj  =  p,  the  covariance  of  z  is  seen  to  be  given  by 


Ef (I  -  E(z]) (I  -  E[z])T] 


Ef(p  -  E[p])(p  -  E[p]Tl 


=  Ef (p  -  p) (p  -  p)  1 


E[[~  E  (P  +  v.)  -  p][i  E  (p  +  v  ) 

n  i  1  n  i=l  1 

*4  C^i) 

n  n  vi=l  '  \j=l  3  ' 

4j  E  E[vt  v*)  ,  since  E(vi  vT] 

n  i=l  J 


P]T] 


=  0  for  i/j 


R__ 

n 

At 

r°  — 2 
R  At 


(52) 


where  R° ,  the  covariance  of  the  measurement  noise  v  ,  is  assumed 
constant  over  the  interval  (tk»  t  +  At], 


In  the  computer  simulations  done  for  this  report,  the  presmoothing 

of  the  measurements  was  simulated  by  scaling  the  measurement  noise  by 

the  factor  \f  At  /At ' .  This  corresponds  to  decreasing  the  measurement 
v  o 

noise  covariance  by  the  factor  AtQ/At.  It  should  b<  noted  that  this 
presmoothing  assumption  is  overly  optimistic. 


The  estimation  error  plots  of  Figs.  6  through  11  illustrate  the 
effects  of  data  presmoothing  upon  the  extended  Kalman  filter’s  performance 
for  Case  1.  Figures  6,  8,  and  10  show  the  effects  of  increasing  the 
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FIG.  6  ESTIMATION  ERRORS  FOR  FULLY-IMPLEMENTED  KALMAN  FILTER  —  CASE 
At  =  0.10  sec  WITHOUT  PRESMOOTHING 
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FIG.  7  ESTIMATION  ERRORS  FOR  FULLY-IMPLEMENTED  KALMAN  FILTER  —  CASE 
At  =  0.10  sec  WITH  PRESMOOTHING 
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FIG.  8  ESTIMATION  ERRORS  FOR  FULLY-IMPLEMENTED  KALMAN  FILTER  —  CASE 
At  0.25  sere  WITHOUT  PRESMOOTHING 
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FIG.  9  ESTIMATION  ERRORS  FOR  FULLY-IMPLEMENTED  KALMAN  FILTER  —  CASE 
At  =  0.25  sec  WITH  PRESMOOTHING 
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FIG.  10  ESTIMATION  ERRORS  FOR  FULLY-IMPLEMENTED  KALMAN  FILTER  —  CASE 
At  =  0.50  sec  WITHOUT  PRESMOOTHING 
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FIG.  11  ESTIMATION  ERRORS  FOR  FULLY-IMPLEMENTED  KALMAN  FILTER  —  CASE 
At  =  0.50  sec  WITH  PRESMOOTHING 


sample  interval  (decreasing  the  data  rate)  from  At  =  0.05  sec  to  At  =  0.10, 
0.25,  and  0.50  sec,  respectively,  when  no  data  presmoothing  is  done.  In 
other  words,  the  filter's  data  rate  is  decreased  by  processing  fewer  mea¬ 
surements  in  any  given  time.  Figures  7,9,  and  11  show  the  effects  of 
increasing  the  sample  interval  of  the  filter  from  At  =  0.05  sec  to  At  =  0.10, 
0.25,  and  0.50  sec,  respectively,  while  also  using  the  increased  time  be¬ 
tween  samples  to  presmooth  the  measurement  data.  The  presmoothing  assump¬ 
tions  described  above  are  used  for  the  results  presented  in  Figs.  7,  9, 
and  11. 

Figures  6,  8,  and  10  indicate  that  if  presmoothing  is  not  done  the 
estimation  errors  become  increasingly  larger  as  the  filter's  sample 
interval  At  is  increased.  However,  Figs.  7,  9,  and  11  indicate  that 
if  the  data  is  presmoothed  the  position  and  velocity  estimation  errors 
generally  decrease  with  increasing  At.  Thus,  the  reduction  of  the 
measurement  noise  through  presmoothing  is  beneficial  in  position  and 
velocity  estimation. 

The  random  disturbance  covariance  Q(k)  given  by  Eq.  (51)  is 
seen  to  increase  linearly  with  increasing  At.  For  ballistic  coefficient 
and  acceleration  estimation,  this  increase  in  Q(k)  apparently  domi¬ 
nates  the  decrease  in  the  measurement  noise  covariance  R(k)  when  the 
data  is  presmoothed.  That  is,  the  errors  in  the  $  and  acceleration 
estimates,  as  shown  in  Figs.  7,  9,  and  11,  increase  with  increasing  At, 
even  with  presmoothing  of  the  data.  Hence,  increasing  the  sample 
interval  generally  leads  to  a  deterioration  in  ballistic  coefficient 
and  acceleration  estimation.  This  is  to  be  expected  since  the  square 
root  of  Q??(k)  [see  Eq.  (51)],  which  is  the  standard  deviation  in  the 
assumed  random  disturbance  noise  affecting  ,  is  of  the  same  order  of 
magnitude  as  the  actual  values  of  x?;  whereas,  the  standard  deviation 
of  the  position  disturbance  noise  is  zero,  and  the  standard  deviation  of 
the  velocity  disturbance  noise  is  very  small  in  relation  to  the  actual 
velocity  values  [see  Eq.  (51)].  Since  the  ballistic  coefficient  8  is 
derived  from  x?  (i.e.,  8  =  p/x7) ,  the  random  disturbance  affecting 

x7  is  of  major  importance  in  estimating  (3.  Consequently,  since  Q?7(k) 
increases  linearly  with  the  sample  interval  At  with  or  without  pre¬ 
smoothing,  it  follows  that  3  estimation  errors  are  expected  to  increase 
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with  increasing  At.  In  turn,  the  errors  in  ballistic  coefficient  estima¬ 
tion  lead  to  errors  in  the  estimate  of  acceleration,  since  the  aerodynamic- 
drag  accelerations  are  in  error. 

For  Case  1,  the  pairs  of  Figs.  6  and  7,  8  and  9,.  10  and  11  illustrate 
the  extended  Kalman  filter  performance,  with  and  without  data  presmoothing, 
for  sample  intervals  of  At  =  0.10,  0.25,  and  0.50  sec,  respectively.  As 
noted  above,  the  presmoothing  assumption  of  Eq.  (52)  is  optimistic;  and 
for  a  given'  At,  the  filter  estimation  errors  for  actual  presmoothing  will 
be  somewhere  between  the  errors  given  in  the  corresponding  pair  of  figures. 
Thus,  for  the  particular  simulations  presented  here,  we  have  bounded  the 
filter  estimation  errors  for  different  values  of  At. 

Figures  12  a,  b,  and  c  show  the  standard  deviations  in  the;  estimates 
of  position,  velocity,  and  3  as  computed  from  the  diagonal  elements  of 
the  calculated  estimation  error  covariance  [Eq.  (29) j ,  Each  of  the  Figs. 

12  a,  b.  and  c  shows,  two  curves:  Curve  (1)  is  the  calculated  standard 
deviation  when  data  is  net  presmoothed,  and  Curve  (ii)  is  the  calculated 
standard  deviation  when  presmoothing  has  been  perfprmed;  these  plots  are 
for  Case  1,  with  a  sample  interval  of  At  =  0.25  sec.  The  standard 
deviation  in  position  error  is  calculated  as  the.  square  root  of  the  sum 
of  the  first  three  diagonal  elements  of  the  filter  covariance  matrix, 

P(k/k) .  That  is,  the  standard  deviation  of  position  error  is  given  by 

standard  deviation  of  velocity  error 
.  The  standard  deviation  for  the 
error  in  0  estimation  is  calculated  using  a  rough  approximation.  The 
ballistic  coefficient  0  is  not  estimated  directly,  but  is  derived  from 
x^.  Since  is  defined  as  =  p/0  (with  p  =  mass  density  of  the 

atmosphere),  it  follows  that  0  =  p/x^.  To  first  order,  we  have 

§  "  P  =  7772  (*7  "  V  * 

(x?) 

At 

so  that  a  rough  approximation  for  the  standard  deviation  of  8  is 


in  which  P  is  the  seventh  diagonal  element  of  P(k/k)  and  o  -  p(h) , 

The  standard  deviation  plots  of  Fig.  12  a,  b,  and  c  provide  a 
valuable  measure  of  the  filter's  performance.  For  instance,  if  the 


P  +  P  +  P 
11  22  33 


Likewise,  the 


is  calculated  as  vP,,  +  P„  +  P„ 
*  44  55  ( 
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FIG.  12  STANDARD  DEVIATIONS  FOR  FULLY-IMPLEMENTED  KALMAN  FILTER  —  CASE 
At  =  0.25  sec:  (i)- WITHOUT  PRESMOOTHING,  (ii)  WITH  PRESMOOTHING 


magnitude  of  the  position  estimation  error  differs  significantly  from 
the  calculated  standard  deviation  of  the  position  error  [as  computed 
from  P(k/k)j  ,,  then  the  filter  is  doing  the  estimation,  task  poorly — : 
perhaps  because  of  inaccurate  modeling  of  the  missile  dynamics  or 
measurement  noise  in  the  filter  equations.  Also ,  when  presmoothing;  of 
the  measurements  leads  to  a  reduction  in  the  magnitude  of  the  estimation 
errors,  the  calculated  standard  deviations  should  also  show  a  reduction; 
Figure  12  shows  that  the  calculated  standard  deviations  [derived  from 
P(k/k) ]  are  in  correspondence  with  the  actual  estimation  error  magnitudes 
(both  with  and  without  presmoothirig) „>  Also,  the  standard  deviations  in 
position  and  velocity  estimation  errors  are  reduced  (as  expected)  when 
data  Is  presmoothed.  Comparison  of  Figs;  8  and  9  with  Fig;  12  reveals 
that  the  actual  estimation  errors  are  less  than  two  standard  deviations 
for  almost  all  of  the  filtering  time. 

Figures  13  a,  b,  and  c  show  the  standard  deviations'  in  position, 
velocity,  and  8  estimation  errors  for  (i)  a  sample;  interval  of  0.05 
sec,  and  (ii)  a  sample  interval  of  0;50  sec;  Curves  (i)  and  (ii)  are  for 
Case  1  with  no  data  presmoothing.  Again,  the  standard  deviation  plots' 
agree  with  the  actual  estimation  error  magnitudes;  For  instance, 
comparison  of  Figs.  4  and  10  with  Fig.  13  shows  that  the  position  esti¬ 
mation  error  magnitude  increases  when  At  is  increased  from  0.05  sec  to 
0.50  sec,  but  this  is  to  be  expected  since ,,  with  no  presmoothing,  the 
calculated  position  standard  deviation  (Fig.  13  a)  also  increases. 

Figures  14  through  19  show  the  estimation  error  magnitudes  when 
the  sample  interval  length  is  increased  from  At  =  0.05  sec  to  At  =  0.10, 
0.25,  and  0.50  sec  for  Case  2.  Estimation  error  magnitudes  are  pre¬ 
sented  for  the  situation  where  there  is  presmoothing  of  measurements  and 
also  where  there  is  no  presmoothing.  The  conclusions  are  similar  to 
those  drawn  for  Case  1  from  Figs.  6  through  11. 

Figure  20  shows  the  effect  of  presmoothing  upon  the  calculated 
standard  deviations  of  the  estimation  errors  for  a  sample  interval  of 
At  =  0.25  sec  in  Case  2.  Again,  the  conclusions  are  the  same  as  those 
drawn  for  Case  1  from  Fig.  12. 
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FIG.  18  ESTIMATION  ERRORS  FOR  FULLY -IMPLEMENTED  KALMAN  FILTER  CASE 
At  =  0.50  sec  WITHOUT  PRESMOOTHING 


0.50  sec  WITH  PRESMOOTHING 


FIG.  20  STANDARD  DEVIATIONS  FOR  FULLY-IMPLEMENTED  KALMAN  FILTER  —  CASE  2, 
At  *  0.25  sec:  (i)  WITHOUT  PRESMOOTHING,  (ii)  WITH  PRESMOOTHiNG 


Figure  21  shows  the  effect  (for  Case  2)  upon  the  calculated 
standard  deviations  in  the  estimation  errors  of  decreasing  the  data  rate. 
.The  plots  show  the  effect  of  decreasing  the  data  rate  from  20  measure¬ 
ments  per  second  to  2  measurements  per  second— when  there  is  no  measure¬ 
ment  presmoothing.  Again,  the  conclusions  are  the  same  as  those  drawn 
for  Case  1  from  Fig.  13. 

A  recent  study6  discusses  the  sensitivity  of  the  Kalman  filter 
performance  to  changes  in  the  sample  interval  At.  In  particular,  the 
effects  of  data  presmoothing  are  discussed  there.  It  was  shown  that  the 
.presmqothing  assumed  here  tends  to  be  overly  optimistic;  since  any 
change  of  the  noise  assumptions  tends  to  cause  the  filter  performance  to 
deteriorate,  from  its  performance  with  the  assumed  presmoothing  [Eq.  (52) J. 

D.  Performance  of  the  Piecewise-Recursive  Kalman  Filter 

The  figures  and  timing  estimates  given  below  indicate  that  the 
piecewise-recursive  Kalman  filter  (described  in  Sec.  IV)  Is  an  extremely 
attractive  alternative  to  the  fully  implemented  (extended)  Kalnian  filter 
for  real-rtime  endoatmospheric  estimation.  Using  this  filtering  method, 
it  is  possible  to  obtain  large  reductions  in  computation  time  with  little 
sacrifice  in  estimation  accuracy.  With  the  Univac  1108  computer,  it  is 
possible  (as  shown  below)  to  filter  measurements  at  a  speed  approximately 
25  times  faster  than  real  time  for  endoatmospheric  estimation,  and  yet 
obtain  accuracy  approaching  that  of  the  fully  implemented  Kalman  filter. 

Since  the  fully  implemented  Kalman  filter  is  approximately  five 
times  faster  than  real-time  (see  Sec.  Ill  C)  on  the  Univac  1108  for 
endoatmospheric  trajectory  estimation  (At  =  0.05  sec),  this  means  that 
the  piecewise-recursive  Kalman  filter  is  about  five  times  faster  than 
the  fully  implemented  Kalman  filter  (the  relative  speed  of  the  two 
filters  is  computer-independent).  With  the  Univac  1108  Computer,  a 
total  time  of  0.010  sec  is  needed  to  process  each  measurement  and  com¬ 
pute  an  estimate  of  the  state  of  the  ballistic  missile  (as  well  as  the 
estimation  error  covariance)  using  the  fully  implemented  Kalman  filter. 
When  the  piecewise-recursive  Kalman  filter  is  used  for  the  estimation, 
two  timing  figures  should  be  noted.  This  arises  from  the  fact  (as  dis¬ 
cussed  in  Sec.  IV)  that  there  are  two  different  modes  of  operation  for  the 
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FIG.  21  STANDARD  DEVIATIONS  FOR  FULLY-IMPLEMENTED  KALMAN  FILTER  —  CASE  2, 
WITHOUT  PRESMOOTHING:  (i)  At  -  0.05  sec,  (ii)  At  .  0.50  sec 


piecewise-recursive  Kalman  filter:  in  the  first  mode,  the  filter  cal¬ 
culates  a  weighting  matrix,  two  error  covariance  matrices,  and  a  state 
estimate;  in  the  second  mode,  the  filter  only  calculates  a  state  esti¬ 
mate  without  the  weighting  matrix  or  error  covariance  computations  being 
performed. 

From  the  simulations  on  the  Univac  1108  for  the  endoatmospheric 
cases,  it  was  found  that  in  the  first  mode,  the  piecewise-recursive 

Kalman  filter  used  an  average  of  0.0173  sec  per  iteration,  while  it  used 

an  average  of  only  0.0007  sec  per  iteration  in  the  second  mode.  Since 

on  the  Univac  11Q8  a  filter  iteration  takes  0.0100  sec  using  the  fully 

implemented  (extended)  Kalman  filter,  this  means  that  in  the  first  mode, 
the  piecewise-recursive  Kalman  filter  is  approximately  1.7  times  slower 
than  the  fully  implemented  Kalman  filter;  however,  in  the  second  mode, 
the  piecewise-recursive  Kalman  filter  is  about  14  times  faster  than  the 
fully  implemented  Kalman  filter.  The  piecewise-recursive  Kalman  filter 
encounters  a  first-mode  calculation  once  every  At  seconds;  and  if 
N  =  At/A1 ,  it  encounters  a  second-mode  calculation  N-l  times  during  a 
time  span  of  At  seconds.  Thus,  for  large  values  of  N,  the  speed  of 
the  piecewise-recursive  Kalman  filter  in  the  second  mode  predominates 
over  that  of  the  first  mode.  For  the  endoatmospheric  trajectories,  the 
average  iteration  time  of  the  piecewise-recursive  Kalman  filter  as  a 
function  of  N  (for  the  values  of  N  encountered  in  the  computer  simu¬ 
lations)  is  given  in  Table  III, 

Using  the  values  in  Table  III,  one  can  calculate  the  average  speed 
of  the  piecewise-recursive  Kalman  filter  relative  to  real  time  for  dif¬ 
ferent  At  and  At.  For  instance,  if  At  =0.5  sec  and  At  =  0.10  sec, 
then  the  speed  of  the  filter  is  about  25  times  faster  than  real  time  (on 
the  Univac  1108);  if  At  =  0.5  sec  and  At  =  0.05  sec,  then  the  speed  of 
the  filter  is  about  20  times  faster  than  real  time  (on  the  Univac  1108). 
Thus,  the  piecewise-recursive  Kalman  filter  can,  on  the  average,  process 
incoming  measurements  20  to  25  times  faster  than  they  are  received. 
Effectively  then,  one  Kalman  filter  of  this  type  could  give  state  esti¬ 
mates  for  20  to  25  endoatmospheric  ballistic  targets  in  real  time. 
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ITERATION  TIMES  FOR  PIECEWISE-RECURSiVE  KALMAN  FILTER 
(ENDOATMOSPHERIC  CASES) 


N  *  — 

At 

AVERAGE  FILTER 
ITERATION  TIME 
(sec) 

ITERATION  TIME 

IN  FIRST  MODE 
(sec) 

ITERATION  TIME 

IN  SECOND  MODE 
(sec) 

3 

0.0061 

- 

C 

o' 

6 

0.0035 

0.0173 

0.0007 

10 

0.0024 

15 

0.0019 

20 

mam 

30 

0.0013 

Now,  let  us  discuss  the  accuracy  of  the  state  estimation  using  the 
piecewise-recursive  Kalman  filter.  It  should  be  pointed  out  that  in  all 
of  the  figures  of  this  section  (for  both  Case  1  and  Case  2)  the  filter 
was  allowed  to  operate  as  a  fully  implemented  Kalman  filter  for  the  first 
four  iterations;  this  operation  serves  as  an  initialization  for  the 
piecewise-recursive  Kalman  filter.  The  estimation  errors  obtained 
through  simulation  for  Case  1  are  exhibited  in  Figs.  22  through  27. 
Figures  22,  23,  and  24  show  the  estimation  errors  for  At  =  0.5,  1.0, 
and  1.5  sec,  respectively,  when  the  sample  interval  is  At  =  0.05  sec. 
Figures  25,  26,  and  27  show  the  estimation  errors  for  At  =  0.5,  1.0, 
and  1.5  sec,  respectively,  when  the  sample  interval  is  At  =  0.10  sec, 
and  there  is  no  data  presmoothing.  The  estimation  errors  exhibited  in 
Figs.  22  through  27  indicate  that,  for  Case  1,  the  piecewise-recursive 
Kalman  filter  has  acceptable  performance  for  At  as  large  as  1.5  sec, 
with  At  =0.05  or  0.10  sec. 


Recall  that  the  error  magnitudes  plotted  in  Fig.  4  for  the  fully 
implemented  (extended)  Kalman  filter  with  a  sample  interval  of  At  =  0.05 
sec  show  that  after  3.5  sec  of  filtering,  the  error  in  estimated  position 
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FIG.  23  ESTIMATION  ERRORS  FOR  PIECEWISE-RECURSIVE  KALMAN  FILTER  —  CASE 
At  =  0.05  sec  AND  At  =  1.0  sec 


5000 


FIG.  24  ESTIMATION  ERRORS  FOR  PIECEWISE-RECURSIVE  KALMAN  FILTER  —  CASE 
At  =  0.05  sec  AND  A r  =  1.5  sec 


Is  less  than  830  ft,  the  error  in  estimated  velocity  is  less  than 

2 

285  ft/sec,  and  the  error  in  estimated  8  is  less  than  665  lb/ft  . 

Figures  22  through  27  indicate  that  for  Case  1  the  piecewise-recursive 

Kalman  filter  estimation  errors  do  not  converge  as  rapidly  as  those  of 

the  fully  implemented  Kalman  filter;  however,  in  the  worst  case,  the 

position  estimation  errors  are  less  than  830  ft  after  6.2  sec,  the 

velocity  estimation  errors  are  less  than  285  ft/sec  after  8.1  sec, 

.  2 

ahd  the  8  estimation  errors  are  less  than  665  lb/ft  after  12.0  sec. 

Figure  28  shows  the  standard  deviations  (for  Case  1)  of  position, 
velocity,  and  8  estimation  errors  [as  calculated  from  P(k/k) ]  for 
(i)  the  fully  implemented  Kalman  filter  with  a  sample  interval, 

At  =  0.10  sec,  and  (ii)  the  piecewise-recursive  Kalman  filter  with 
At  =  0.10  sec  and  At  =  1.0  sec.  These  two  curves  indicate  that  the 
assumption  made  in  Eqs.  (36)  and  (37)  of  Sec.  IV  concerning  the  approxi¬ 
mation  of  the  filtering  dene  by  the  piecewise-recursive  Kalman  filter 
is  valid.  Specifically,  we  had  considered  the  filtering  done  with  a 
fixed  weighting  matrix  as  being  equivalent  to  measurement  presmoothing, 
so  that  when  the  estimation  error  covariance  matrix  was  updated  [see 
Eq.  (37)],  the  calculation  used  a  measurement -noise  covariance  reduced 
to  R(k) At/AT.  The  validity  of  this  assumption  is  verified  in  Fig,  28, 
since  each  new  calculation  of  the  estimation  error  covariance  at  inter¬ 
vals  of  At  seconds  gives  a  decrease  to  nearly  the  value  of  the 
covariance  matrix  obtained  with  the  use  of  the  fully  implemented  Kalman 
filter  when  it  is  calculated  at  intervals  of  At  seconds. 

Figures  29  through  34  show  the  estimation  error  magnitudes  obtained 
through  computer  simulation  for  Case  2.  Figures  29,  30,  and  31  show  the 
estimation  errors  for  At  =  0.5,  1.0,  and  1.5  sec,  respectively,  when 
the  sample  interval  At  =  0,05  sec.  Figures  32,  33,  and  34  show  the 
estimation  errors  for  At  =  0.5,  1.0,  and  1.5  sec,  respectively,  when 
the  sample  interval  is  At  =  0.10  sec.  For  Case  2,  we  see  from  Figs.  29 
through  34  that  for  either  At  =  0.05  sec  or  At  =  0.10  sec,  the  filter 
performance  deteriorates  as  At  increases  from  0.5  sec  to  1.5  sec. 


FIG.  28  STANDARD  DEVIATIONS  FOR  PIECEWISE-RECURSIVE  KALMAN  FILTER  —  CASE 


FIG.  31  ESTIMATION  ERRORS  FOR  PIEC.EWISE-RECURSIVE  KALMAN  FILTER  ~  CASE  2, 
At  =  0.05  sec  AND  At  =  1.5  sec 


Figure  35  shows  the  standard  deviations  (for  Case  2)  of  position, 
velocity,  and  8  estimation  errors  [as  calculated  from  P(k/k)]  for 
(i)  the  fully  implemented  Kalman  filter  with  a  sample  interval  of 
At  =  0.10  sec,  and  (ii)  the  piecewise-recursive  Kalman  filter  with 
At  =  0.10  sec  and  At  -  1.0  sec.  As  in  Case  1,  Fig.  35  indicates  that 
the  assumption  concerning  the  approximation  of  the  measurement-noise 
covariance  in  the  estimation  error  covariance  for  the  piecewise- 
recursive  filter  is  valid  for  Case  2. 

In  summary,  the  computer  simulations  indicate  that  for  endoatmospheric 
estimation  the  accuracy  of  the  piecewise-recursive  Kalman  filter  is 
comparable  with  that  of  the  fully  implemented  Kalman  filter,  and  that  it 
offers  significant  savings  in  computation  time. 


VI  NUMERICAL  RESULTS  FOR  EXOATMOSPHERIC  ESTIMATION 

A.  Description  of  the  Test  Case  and  Filter.  Design  Parameters 

Data  from  a  representative  exoatmospheric  trajectory  generated  by 
simulation  using  an  accurate  model  of  reentry  dynamics  was  used  to 
evaluate  the  performance  of  the  piecewise-recursive  Kalman  filter 
described  in  Sec.  IV.  The  computer  program  used  in  generating  this 
trajectory  was  developed  by  M.I.T.  Lincoln  Laboratory. 

In  this  report,  the  exoatmospheric  test  trajectory  will  be  referred 
to  as  Case  3.  It  is  a  minimum-energy  ellipse  having  a  range  of  approxi- 

7 

mately  2.6  X  10  ft  (4300  nmi)  from  launch  to  impact,  with  a  reentry 
angle  of  about  23°.  This  trajectory  is  estimated  from  a  range  of 

e  g 

9.4  X  10  ft  (1600  nmi)  to  a  range  of  0,4  X  10  ft  (66  nmi),  which 
corresponds  to  an  altitude  of  approximately  370,000  ft  for  the  particular 
geometry  of  this  case. 

At  altitudes  greater  than  300,000  ft,  the  atmospheric  drag  is 
negligible;  hence,  in  the  filter's  model  for  the  equations  of  the  bal¬ 
listic  missile,  the  drag  terms  have  been  dropped  [see  Eqs.  (1)  in  Sec. 

II  A],  Thus,  the  state  variable  x?,  which  contains  8,  can  also  be 
eliminated;  of  course,  the  dimensions  of  the  system  matrices  ($ ,  H,  W, 
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FIG.  35  STANDARD  DEVIATIONS  FOR  PIECEWISE-RECURSIVE  KALMAN  FILTER  —  CASE  2, 


Q,  and  P)  are  modified  accordingly.  Due  to  numerical  difficulties  that 
were  encountered  in  calculating  the  covariance  matrices,  the  filter  was 
programmed  in  double  precision  for  exoatmospheric  estimation. 

The  filtering  computations  start  with  an  initial  state  estimate 
£(°/0) »  which  is  now  a  six-dimensional  vector,  and  its  error  covariance 
P(0/0) ,  The  initial  state  estimates  for  the  exoatmospheric  case  are 
given  in  Table  IV. 


Table  IV 

INITIAL  CONDITIONS  FOR  EXOATMOSPHERIC  TRAJECTORY 


X1 

(ft  x  106) 

x2 

(ftx  106) 

x3 

(ft  x  106) 

X4 

(ft/ sec) 

x5 

(ft/sec) 

x6 

(ft/sec) 

CASE  3 

ACTUAL:  x  (0) 

-6.01 

7.04 

1.73 

11500 

-12200 

3020 

ESTIMATE:  x  (0/0) 

-6.10 

7.14 

1.63 

12300 

-11420 

3268 

{ERROR] 

0.09 

0.10 

0.10 

800 

780 

248 

The  initial  error  covariance  (for  the  exoatmospheric  trajectory) 
was  determined  experimentally  to  be  the  6x6  diagonal  matrix 


10 


12 


P(0/0)  = 


10 


12 


10' 


(53) 


The  measurement-noise  covariance  R(k)  is  based  upon  a  model  of 
actual  radar-noise  characteristics  of  the  PAR;  the  model  was  obtained 
by  making  suitable  modifications  of  the  results  presented  in  Ref.  9. 
This  radar  model  is  used  to  generate  the  measurement  noise  of  Eq.  (12) 
in  the  simulation.  Observations  may  be  taken  as  often  as  every  second 
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for  the  exoatmospheric  case,  but  the  filter  may  process  fewer  measure¬ 
ments  than  are  taken  (i.e.,  At  £  1.0  sec). 

The  random  disturbance  covariance  Q(k)  [see  Eqs.  (15a)  and  (15b)] 
compensates  for  model  inaccuracies  as  described  in  Sec.  II  A.  This 
covariance  was  determined  experimentally  to  be  the  6x6  diagonal 
matrix 


B.  Performance  of  the  Fully  Implemented  Extended  Kalman  Filter 

For  the  exoatmospheric  trajectory,  the  standard  used  to  evaluate 
the  performance  of  the  piecewise-recursive  Kalman  filter  is  the  fully 
implemented  (extended)  Kalman  filter  with  a  sample  interval  of  At  =  1.0 
sec.  The  recursive  equations  for  the  extended  Kalman  filter  are  given 
in  Sec.  Ill  A,  with  the  system  matrices  (§ ,  H,  V/,  Q,  and  P)  modified 
appropriately.  This  filter  was  used  to  estimate  the  ballistic  trajectory 
for  Case  3,  using  the  simulated  trajectory  and  measurement  noise.  The 
filter  yields  state  estimates  x(k/k)  at  sample  intervals  of  At  =  1.0 
sec. 


0  100  200  300  400  500 

TIME — sec 


TC*M88«800 

F!G.  36  MEASUREMENT  NOISE  —  CASE  3 
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The  magnitude  of  the  measurement  error  in  position  is  shown  in  Fig.  36 
for  Case  3.  In  order  to  relate  range  and  time,  the  measurement  error  is 
plotted  as  a  function  of  both  variables;  in  all  other  figures,  the  plots 
are  a  function  of  time  only.  Figures  37  a,  b,  and  c  show  the  magnitude 
of  the  estimation  errors  for  position,  velocity,  and  acceleration, 
respectively,  for  Case  3.  It  should  be  noted  that  acceleration  is  not 
computed  directly  by  the  filter  but  is  calculated  from  the  state  estimate 
using  Eqs.  (1)  of  Sec.  II  A.  Additionally,  as  noted  above,  3  is  not 
estimated  for  the  exoatmospheric  case. 

From  Fig.  37,  it  can  be  seen  that  after  30  sec  of  filtering,  the 
magnitude  of  the  error  in  estimated  position  is  less  th^n  6,000  ft  (for 
Case  3,  the  initial  error  in  position  is  167,000  ft);  after  170  sec,  the 
position  error  is  less  than  2,000  ft.  After  30  sec  of  filtering  the 
magnitude  of  the  error  in  estimated  velocity  is  less  than  200  ft/sec 
(for  Case  3,  the  initial  error  in  velocity  is  almost  1200  ft/sec);  after 
100  sec,  the  velocity  error  is  less  than  50  ft/sec. 

C.  Performance  of  the  Fully  Implemented  Extended  Kalman  Filter  with 

Different  Data  Rates 

As  discussed  in  Sec.  V  C,  the  simplest  scheme  for  improving  the 
real-time  capability  of  the  Kalman  filter  is  simply  to  decrease  the  data 
rate.  If  the  sample  intervals  are  lengthened  and  the  additional  time 
between  samples  is  not  used  to  presmooth  the  measurements  (i.e.,  the 
intermediate  measurements  are  discarded)  ,  then  obviously  the  filter 
performance  will  deteriorate  to  some  extent.  The  present  study  also 
considers  the  situation  where  there  is  presmoothing  of  the  measurements. 

In  practice,  data  is  presmoothed  by  curve-fitting  with  polynomials  that 
are  not  precise  representations  of  missile  trajectories. 

The  assumption  of  data  presmoothing  implies  that  in  the  filter 
equations  the  measurement -noise  covariance  should  be  reduced  accordingly. 

In  the  computer  simulations  done  for  this  report,  the  presmoothing  of 
the  measurements  was  simulated  by  scaling  the  measurement  noise  by  the 
factor  J  AtQ/At * .  This  corresponds  to  decreasing  the  measurement-noise 
covariance  by  the  factor  AtQ/At.  A  justification  for  the  approximation 
is  posented  in  Sec.  V  C,  although  it  should  be  noted  that  this  pre¬ 
smoothing  assumption  is  overly  optimistic. 
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FIG.  37  ESTIMATION  ERRORS  FOR  FULLY-IMPLEMENTED  KALMAN  FILTER  —  CASE  3,  At  =  1.0  sec 
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2.0  sec  WITHOUT  PRESMOOTHING 
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FIG.  39  ESTIMATION  ERRORS  FOR  FULLY-IMPLEMENTED  KALMAN  FILTER  —  CASE  3, 
At  •=•  2.0  sec  WITH  PRESMOOTHING 
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5.0  sec  WITHOUT  PRESMOOTHING 


The  estimation  error  plots  of  Figs.  38  through  41  illustrate  the 
effects,  of  data  presmoothing  upon  the  extended  Kalman  filter's  perfor-r 
mance  for  Case  3..  Figures  39  and  41  show,  the  effects  of  increasing  the 
sample  interval  (decreasing  the  data  rate)  from  At  =  1.0  sec  to 
At  =2.0  and  5.0  sec,  respectively,,  while  also  using  the  increased  time 
between  samples  to  presmooth  the  measurement  data.  The  presmoothing 
assumptions  described  above  are  used  for  the  results  presented  in 
Figs.  39  arid  41. 

Figures  38  and  40  indicate  that  without  pre smoothing,  the  estimation 
errors  become  increasingly  larger  as  the  filter's  sample  interval  At 
is  increased.  However,  Figs.  39  and  41  indicate  that  if  the  data  is 
presmoothed  the  position  and  velocity  estimation  erroi-s  generally 
decrease  with  increasing  At>  Thus,  the  reduction  of  the  measurement 
noise  through  presmoothing  is  beneficial  in  position  and  velocity  esti¬ 
mation. 

The  pairs  of  Figs.  38  and  39,  40  and  41  illustrate  the  extended 
Kalman-f liter  performance,  with  arid  without  data  presmoothing ,  for 
sample  intervals  of  At  =  2.0  arid  5 ;0  sec,  respectively.  As  rioted 
above,  the  presmoothing  assumption  is  optimistic  arid  for  a  given  At, 
the  filter  estimation  errors  for  actual  presmoothing  will  be  somewhere 
between  the  errors  given  in  the  corresponding  pair  of  figures.  Thus, 
for  the  particular  simulations  presented  here,  we  have  bounded  the 
filter  estimation  errors  for  different  values  of  At. 

Figures  42  a  and  b  show  the  standard  deviations  in  the  estimates  of 
position  and  velocity  as  computed  from  the  diagonal  elements  of  the 
calculated  estimation  error  covariance  [Eq,  (29)].  Each  of  the  Figs.  42 
a  and  b  shows  two  curves:  Curve  (i)  is  the  calculated  standard  devia¬ 
tion  when  there  is  no  data  presmoothing,  and  Curve  (ii)  is  the  calculated 
standard  deviation  when  there  is  presmoothing;  these  plots  are  for 
Case  3,  with  a  sample  interval  of  At  =5.0  sec.  The  standard  deviation 
in  position  error  is  calculated  as  the  square  root  of  the  sum  of  the 
first  three  diagonal  elements  of  the  filter  covariance  matrix,  P(k/k) . 
That  is,  the  standard  deviation  of  position  error  is  given  by 
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deviation  of  velocity  error 

Figure  42  shows  that  the  calculated  standard  deviations,  [derived 
from  P(k/k)]  are  in  correspondence  with  the  actual  estimation  error 
magnitudes  (both  with  and  without  pre smoothing) .  Also,  the  standard 
deviations  in  position  and  velocity  estimation  errors  are  reduced  (as 
expected)  when  there  is  presmoOthing. 

Figures  43  a  and  b  show  the  standard  deviations  in  position  and 
velocity  estimation  errors  for  (i)  a  sample  interval  of  1.0  sec,  and 
(ii)  a  sample  interval  of  5.0  sec;  Curves  (i)  and  (ii)  are  for  Case  3 
with  no  data  presmoothing.  Again,  the  standard  deviation  plots  agree 
with  the  actual  estimation  error  magnitudes.  For  instance*  comparison 
of  Figp  .  37  and  40  with  Fig.  43  shows  that  the  position  estimation  error 
magnitude  increases  when  At  is  increased  from  1.0  -sec  to  5.0  sec,  but 
this  is  to  be  expected  since,  with  no  presmoothing,  the  calculated 
position  standard  deviation  [Fig;  43  a]  also  increases. 

D.  Performance  of  the  PiecewiserRecursive  Kalman  Filter 

The  figures  and  timing  estimates  given  below  indicate  that  the 
piecewise-recursive  Kalman  filter  (described  in  Sec,  IV)  is  an  extremely 
attractive  alternative  to  the  fully  implemented  (extended)  Kalman  filter 
for  real-time  exoatmospheric  estimation.  Using  this  filtering  method, 
it  is  possible  to  obtain  large  reductions  in  computation  time  with  little 
sacrifice  in  estimation  accuracy.  With  the  Univac  1108  computer,  it  is 
possible  (as  shown  below)  to  filter  measurements  at  a  speed  approximately 
500  times  faster  than  real  time  for  exoatmospheric  estimation,  and  yet 
obtain  accuracy  approaching  that  of  the  fully  implemented  Kalman  filter. 

Since  the  fully  implemented  Kalman  filter  is  approximately  100  times 
faster  than  real-time  (see  Sec.  Ill  C)  on  the  Univac  1108  for  exoatmos¬ 
pheric  trajectory  estimation  (At  =  1.0  sec) ,  this  means  that  the  piecewise- 
recursive  Kalman  filter  is  about  5  times  faster  than  the  fully  implemented 
Kalman  filter  (the  relative  speed  of  the  two  filters  is  essentially 
computer-independent).  With  the  Univac  1108  Computer,  a  total  time  of 
0.011  sec  is  needed  to  process  each  measurement  and  compute  an  estimate 


J  P.  +  P,  +  P__  .  Likewise,  the  standard 

V  ii  12  1  ■■■*  '■  ■■■■  | 

is  calculated  as  J  P ...  +  P**.  +  P-  . 

V  44  55  66 
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of  the  state  of  the  ballistic  missile  (as  well  as  the  estimation  error 
covariance)  using  the  fully  implemented  Kalman  filter.  When  the 
piecewlse-recursive  Kalman  filter  is  used  for  the  estimation,  two 
timing  figures  should  be  noted.  This  arises  from  the  fact  (as  discussed 
in  Sec.  IV)  that  there  are  two  different  modes  of  operation  for  the 
piecewise-recursive  Kalman  filter:  in  the  first  mode,  the  filter  cal¬ 
culates  a  weighting  matrix,  two  error  covariance  matrices,  and  state 
estimate;  in  the  second  mode,  the  filter  calculates  only  the  state 
estimate,  without  the  weighting  matrix  or  error  covariance  computations 
being  performed. 

From  the  simulations  on  the  Univac  1108  for  the  exoat mospheric 
case,  it  was  found  that  in  the  first  mode,  the  piecewise-recursive 
Kalman  filter  used  an  average  of  0.0200  sec  per  iteration,  while  it 
used  an  average  of  only  0.0010  sec  per  iteration  in  the  second  mode. 

Since  a  filter  iteration  on  the  Univac  1108  takes  0.0110  sec  using  the 
fully  implemented  (extended)  Kalman  filter,  this  means  that,  in  the 
first  mode,  the  piecewise-recursive  Kalman  filter  is  approximately  1.8 
times  slower  than  the  fully-implemented  Kalman  filter;  however,  in  the 
second  mode,  the  piecewise-recursive  Kalman  filter  is  about  11  times 
faster  than  the  fully  implemented  Kalman  filter.  The  piecewise-recursive 
Kalman  filter  encounters  a  first-mode  calculation  once  every  At  seconds; 
and  if  N  »  At/A1 ,  it  encounters  a  second- lode  calculation  N-l  times 
during  a  time  span  of  At  seconds.  Thus,  for  large  values  of  N,  the 
speed  of  the  piecewise-recursive  Kalman  filter  in  the  second  mode  pre¬ 
dominates  over  that  of  the  first  mode.  For  the  exoatmospheric  trajectory, 
the  average  iteration  time  of  the  piecewise-recursive  Kalman  filter  as  a 
function  of  N  (for  the  values  of  N  encountered  in  the  computer 
simulations)  is  given  in  Table  V. 

Using  the  values  in  Table  V,  one  can  calculate  the  average  speed  of 
the  piecewise-recursive  Kalman  filter  relative  to  real  time  for  different 
At  and  At.  For  instance,  if  AT  =  20  sec,  and  At  =  1.0  sec,  then  the 
speed  of  the  filter  is  about  500  times  faster  than  real  time  (on  the 
Univac  1108).  Thus,  with  =  20  sec  and  At  =  1.0  sec,  the  piecewise- 
recursive  Kalman  filter  can,  on  the  average,  process  incoming  measurements 
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Table  V 


ITERATION  TIMES  FOR  PIECEWISE-RECURSIVE  KALMAN  FILTER 
(EXOATMOSPHERIC  CASES) 


N  -  At 

AVERAGE  FILTER 
ITERATION  TIME 
(sec) 

ITERATION  TIME 

IN  FIRST  MODE 
(sec) 

ITERATION  TIME 

IN  SECOND  MODE 
(sec) 

5 

0.0048 

10 

0.0029 

20 

0.0020 

0.0200 

0.0010 

50 

0.0014 

100 

0.0012 

500  times  faster  than  they  are  received.  Effectively  then,  one  Kalman 
filter  of  this  type  could  give  state  estimates  for  500  exoatmospheric 
ballistic  targets  in  real  time. 

Now,  let  us  discuss  the  accuracy  of  the  state  estimation  using  the 
piecewise-recursive  Kalman  filter.  It  should  be  pointed  out  that  in  all 
of  the  figures  of  this  section  (for  Case  3)  the  filter  was  allowed  to 
operate  as  a  fully  implemented  Kalman  filter  for  the  first  10  iterations; 
this  operation  serves  as  an  initialization  for  the  piecewise-recursive 
Kalman  filter.  The  estimation  errors  obtained  through  simulation  for 
Case  3  are  exhibited  in  Figs.  44  through  49.  Figures  44,  45,  and  46 
show  the  estimation  errors  for  At  =  20,  50,  and  100  sec,  respectively, 
when  the  sample  interval  is  At  =  2.0  sec  and  there  is  no  data  presmoothing. 
The  estimation  errors  exhibited  in  Figs.  44  through  49  indicate  that,  for 
Case  3,  the  piecewise-recursive  Kalman  filter  has  acceptable  performance 
for  At  as  large  as  100  sec,  with  At  =  1.0  or  2.0  sec. 

Recall  that  the  error  magnitudes  plotted  in  Fig.  37  for  the  fully 
implemented  (extended)  Kalman  filter  with  a  sample  interval  of  At  =  1.0 
sec  show  that  after  30  sec  of  filtering,  the  error  in  estimated  position 
is  less  than  6,000  ft  and  the  error  in  estimated  velocity  is  less  than 
200  ft/sec.  Figures  44  through  49  indicate  that  for  Case  3,  the 
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FIG.  44  ESTIMATION  ERRORS  FOR  PIECEWISE-RECURSIVE  KALMAN  FILTER  -  CASE  3, 
At  =  1 .0  sec  AND  Ar  =  20  sec 
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FIG.  46  ESTIMATION  ERRORS  FOR  PIECEWISE-RECURS1VE  KALMAN  FILTER  —  CASE 
At  =  1.0  sec  AND  At  --  100  sec 
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FIG.  47  ESTIMATION  ERRORS  FOR  PIECEWISE-RECURSIVE  KALMAN  FILTER  CASE 
At  =2.0  sec  AND  At  =  20  sec 
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2.0  sec  AND  At  =■  50  sec 
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FIG.  49  ESTIMATION  ERRORS  FOR  P1ECEWISE-RECURSIVE  KALMAN  FILTER  —  CASE  3, 
\t  2.0  sec  AND  \r  .  100  sec 


piecewise-recursive  Kalman  filter  estimation  errors  do  not  converge  as 
rapidly  as  those  of  the  fully  implemented  Kalman  filter  ;  however,  in 
the  worst  case  (At  =  1.0  or  2.0  sec  and  At  =  100  sec)  the  position 
estimation  errors  are  less  than  6,000  ft  after  230  sec  and  the  velocity 
estimation  errors  are  less  than  200  ft/sec  (except  for  infrequent 
excursions)  after  200  sec. 

Figure  50  shows  the  standard  deviations  (for  Case  3)  of  position 
and  velocity  estimation  errors  [as  calculated  from  P(k/k) ]  for  (i)  the 
fully  implemented  Kalman  filter  with  a  sample  interval,  At  =2.0  sec, 
and  (ii)  the  piecewise-recursive  Kalman  filter  with  At  =  2.0  sec  and 
At  =  20  sec.  These  two  curves  indicate  that  the  assumption  made  in 
Eqs.  (36)  and  (37)  of  Sec.  IV  concerning  the  approximation  of  the  fil¬ 
tering  done  by  the  piecewise-recursive  Kalman  filter  is  valid.  Specifi¬ 
cally,  we  had  considered  the  filtering  done  with  a  fixed  weighting 
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matrix  as  being  equivalent  to  measurement  presmoothing,  so  that  when  the 
estimation  error  covariance  matrix  was  updated  [see  Eq.  (37)],  the  cal¬ 
culation  used  a  measurement  noise  covariance  reduced  to  R(k)  At/AT. 

The  validity  of  this  assumption  is  verified  in  Fig.  50,  since  each  new 
calculation  of  the  estimation  error  covariance  at  intervals  At  sec 
gives  a  decrease  to  nearly  the  value  of  the  covariance  matrix  obtained 
with  the  use  of  the  fully  implemented  Kalman  filter  when  it  is  calcu¬ 
lated  at  intervals  of  At  sec. 

In  summary,  the  computer  simulations  indicate  that  for  exoatmospheric 
estimation  the  accuracy  of  the  piecewise-recursive  Kalman  filter  is  com¬ 
parable  with  that  of  the  fully  implemented  Kalman  filter,  and  that  it 
offers  significant  savings  in  computation  time. 

VII  CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FUTURE  STUDY 

This  study  has  addressed  itself  to  the  problem  of  real-time  imple¬ 
mentation  of  the  Kalman  filter  for  estimating  ballistic  trajectories. 

As  shown  in  this  study  and  Ref.  1,  the  Kalman  filter  is  an  extremely 
effective  algorithm  for  estimation  of  ballistic  trajectories,  although 
the  computational  requirements  of  the  fully  implemented  Kalman  filter 
are  quite  severe.  On  the  Univac  1108  Computer,  the  fully  implemented 
Kalman  filter  runs  about  five  times  faster  than  real  time  for  the  endo- 
atmospheric  cases  and  about  100  times  faster  than  real  time  for  the 
exoatmospheric  cases. 

In  this  report,  several  approaches  that  may  be  used  to  modify  the 
Kalman  filtering  algorithm  in  order  to  reduce  the  computational  require¬ 
ments  are  described.  The  most  promising  approach  of  those  considered 
is  the  piecewise-recursive  Kalman  filter.,  which  is  described  in  Sec.  IV. 

As  shown  by  the  numerical  results  obtained  from  extensive  computer 
simulations  (see  Secs.  V  and  VI),  the  piecewise-recursive  Kalman  filter 
can  process  measurements  of  a  single  target  at  a  computational  speed 
(on  the  Univac  1108)  that  is  about  20  to  25  times  faster  than  real  time  for 
the  endoatmospheric  cases  and  about  500  times  faster  than  real  time  for 


91 


the  exoatmospheric  cases,  and  yet  obtain  estimation  accuracy  approaching 
that  of  the  fully  implemented  Kalman  filter.  This  increased  filter 
capability  is  invaluable  for  the  real-time  estimation  of  multiple  targets. 

Other  investigations  that  should  be  performed  in  order  to  further 
reduce  the  computation  time  of  the  piecewise-recursive  Kalman  filter  are: 

(1)  Use  a  piecewise-linear  weighting  matrix  in  the  filter, 
rather  than  a  piecewise-constant  weighting  matrix  as 
described  in  Sec.  IV.  This  approach  can  be  readily  in¬ 
corporated  into  the  existing  computer  program  and  should 
yield  significant  reductions  in  computation  time  for 
comparable  filter  performance  by  allowing  larger  values 
of  At. 

(2)  Vary  At  and/or  At  over  different  portions  of  the 
estimation  interval.  This  approach  can  be  placed  upon  a 
rigorous  basis  by  applying  the  sensitivity  results 
obtained  in  Ref.  6,  and  it  also  should  reduce  the  filter's 
average  computation  time  per  iteration. 

(3)  Simplify  the  model  of  the  system  equations  that  is  used 
in  the  filter's  covariance  calculations.  As  shown  in 
Table  I  of  Sec,  III  C,  these  covariance  calculations 
comprise  the  major  portion  of  the  computation  time  per 
filter  iteration. 
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